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SECTION  1 
INTRODUCTION 


This  report  describes  the  development  of  theoretical  and 
numerical  methods  for  use  in  simulating  the  dynamic  response  of 
structures  to  hydrodynamic  loading.  The  particular  problem  of 
interest  is  that  of  soft— body  impacts  on  aircraft  transparencies/ 
which  involves  several  characteristics  which  make  effective 
analysis  guite  difficult.  Among  these  characteristics  are 

-  nonlinear  structural  response  (large  displacements  and 

nonlinear  material  behavior); 

-  structural  failures  which  often  occur  at  times  far  in 

excess  of  the  duration  of  the  impact; 

-  very  large  deformation  of  the  impacting  body; 

-  strong  interaction  between  structural  response  and  the 

magnitude  and  distribution  of  impact  pressures;  and 

-  free-surface  boundary  conditions  on  the  impacting  body. 

The  analytical  techniques  presented  here  address  each  of  these 
problem  areas  to  provide  effective  numerical  solutions  to 
practical  problems  involving  soft  body  impact. 

The  remainder  of  this  Section  describes  the  technical 
aspects  of  the  soft-body  impact  problem  in  more  detail/  and 
outlines  the  general  approach  adopted  in  its  solution. 

Subsequent  Sections  are  devoted  to  presenting  the  theoretical 
development  and  essential  information  about  the  computer 
implementation  of  the  methodology. 

1.1  PROBLEM  OVERVIEW 

Contact  and  impact  situations  are  among  the  most  important 
and  difficult  problems  in  computational  mechanics.  Contact  or 
impact  between  solid  bodies  occurs  in  forming  operations/  wheel- 
rail  contact/  ballistic  impacts/  and  numerous  other  practical 
circumstances.  Fluid-structure  interaction/  contact  between 
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solid  and  fluid  media,  is  of  interest  in  submerged  structures, 
reactor  vessels,  and  turbomachinery.  For  each  of  these  classes 
of  analytical  problems,  solution  procedures  of  many  types  have 
been  developed  and  used  successfully  [1-5].  For  the  case  of 
solid  to  solid  contact,  very  general  analysis  techniques  exist, 
based  upon  both  explicit  and  implicit  integration  methods.  In 
fluid-structure  interaction  problems,  the  solution  procedure 
necessarily  reflects  the  specific  problem  geometry  or  response 
characteristics  rather  strongly. 

Soft-body  impacts,  such  as  bird  impacts  on  aircraft 
transparencies,  involve  the  interaction  of  a  small  volume  of 
solid  material  which  behaves  hydrodynamically  (similar  to  a 
fluid)  with  a  more  massive  structure  whose  response  is  generally 
less  severe.  When  Eulerian  methods  of  description  (based  on  a 
mesh  fixed  in  space)  are  used  for  the  hydrodynamic  solution, 
tracking  of  boundaries  and  free  surfaces  is  difficult. 

Lagrangian  mesh  techniques,  in  which  computational  reference 
points  move  with  the  material,  fail  due  to  the  severe  mesh 
distortion  which  may  occur  within  the  soft  body.  Explicit,  large 
deformation  codes  for  solid  mechanics  applications  [6-8] 
sometimes  provide  the  type  of  analytical  capabilities  which  are 
necessary  to  predict  the  resulting  behavior;  however,  the  need  to 
compute  structural  response  over  relatively  long  times  (several 
milliseconds)  makes  the  use  of  these  techniques  prohibitive  in 
terms  of  cost.  Due  to  the  lack  of  solution  methodology  which  is 
truly  tailored  to  the  soft— body  impact  problem,  ad-hoc  techniques 
abound  [9—11],  While  such  methods  are  economical  and  provide 
useful  data  for  preliminary  design,  their  reliability  leaves  much 
to  be  desired,  and  the  objective  of  replacing  experiment  by 
analysis  remains  unfulfilled. 
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1.2  GENERAL  APPROACH 


The  solution  techniques  presented  here  are  based  upon  the 
finite  element  method  [12],  Due  to  the  widely  differing 
characteristics  of  the  soft  body  and  the  structure,  the  coupled 
impact  problem  is  solved  separately  (but  in  tandem)  for  each 
region,  with  appropriate  interface  conditions  being  imposed  at 
each  stage  of  the  solution. 

The  structural  model  experiences  relatively  small 
deformations,  which  must  be  computed  for  extended  time  periods  in 
some  cases.  For  the  most  part,  structural  response  (as  opposed 
to  wave  propagation)  is  of  primary  interest.  Therefore,  the 
finite  element  solution  for  the  structure  motion  uses 
higher-order  elements  in  conjunction  with  highly  stable,  implicit 
integration  techniques.  This  combination  limits  the  model  to 
reasonable  size  and  permits  the  use  of  time  steps  in  the 
millisecond  range.  In  return  for  this  convenience,  the  implicit 
method  requires  an  iterative  solution  at  each  time  step,  and 
formation  of  a  global  stiffness  matrix.  This  approach  to  the 
dynamic  response  solution  is  common  in  general-purpose  finite 
element  packages,  such  as  ABAQUS  [13],  ADINA  [14],  ANSYS  [15], 
MAGNA  [16]  and  MARC  [17]. 

The  behavior  of  the  soft  body  region  is  quite  different  in 
character  from  the  structural  response.  Here  the  deformations 
are  large  and,  due  to  the  low  material  stiffness,  depend  upon 
wave  propagation  effects  as  well  as  overall  deformation  modes. 

The  response  is  essentially  hydrodynamic;  that  is,  pressure 
levels  are  typically  many  times  greater  than  the  strength  of  the 
material.  For  this  region,  the  numerical  solution  is  based  upon 
explicit  techniques  which  permit  a  much  more  detailed  treatment 
of  the  problem  physics.  Due  to  the  free  surface  which  makes  up 
most  of  the  boundary  of  the  body,  a  Lagrangian  mesh,  which  moves 
with  the  material,  is  used.  The  severe  mesh  distortion  which 
occurs  in  some  situations  is  resolved  by  occasional  "rezoning"; 
that  is,  complete  regeneration  of  the  mesh.  While  mesh  rezoning 
is  often  a  source  of  difficulty  in  large  deformation  analyses, 
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the  present  technique  includes  automatic  rezoning  operations 
which  are  invoked  when  preset  tolerances  are  exceeded,  and  do  not 
require  stopping  and  restarting  the  solution. 

Appropriate  interface  conditions  between  the  structure  and 
soft  body  meshes  consist  of  kinematic  conditions  (compatibility 
of  displacements  and  velocities)  and  force/impulse  conditions 
(interface  pressures  and  impulses).  In  the  structural  model, 
forces  exerted  by  the  hydrodynamic  mesh  are  applied  directly  to 
compute  incremental  displacements;  for  the  soft  body,  imposed 
displacement  conditions  are  provided  by  the  structure  model,  and 
interface  forces  develop  as  natural  boundary  conditions.  For  the 
c on ta c t/i mpa c t  calculations,  the  algorithm  presented  in  Reference 
[18]  is  used.  Since  stable  time  steps  for  the  hydrodynamic  mesh 
are  much  smaller  than  for  the  structure  mesh,  a  partitioning 
method  [19]  is  needed  to  enforce  the  proper  flow  of  information 
between  the  two  distinct  portions  of  the  mesh. 

The  primary  development  areas  in  the  work  described  are  the 
hydrodynamic  solution  and  the  interface  techniques.  These  are 
suitable  for  integration  into  most  implicit,  nonlinear  structural 
analysis  programs.  In  the  present  effort,  we  have  chosen  one 
particular  program,  MAGNA  (Materially  And  Geometrically  Nonlinear 
Analysis)  in  which  to  implement  the  coupled  impact  solution. 
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1.3  NOTATION  CONVENTIONS 

The  notation  used  in  the  body  of  this  report  conforms  to  the 
following  rules  and  conventions: 

-  Einstein's  summation  convention  applies  unless  otherwise 
noted;  that  is,  repeated  indices  imply  summation  over  the  range 
of  the  index; 

-  lower  case  subscripts  are  used  for  reference  to  spatial 
dimensions  and  have  a  range  of  three; 

-  upper  case  subscripts  refer  to  the  nodes  of  a  finite 
element,  and  therefore  have  a  range  of  N  for  an  N-noded  element; 

-  material  time  derivatives  are  indicated  by  a  superimposed 
dot,  C); 

-  a  comma  represents  partial  differentiation  with  respect 
to  the  spatial  coordinate  which  follows  (e.g.,  u  .  -  9u/3x.); 
where  both  time  and  spatial  derivatives  occur,  as  in  u  .,  the 

f  1 

material  time  derivative  is  evaluated  first. 
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SECTION  2 

THEORETICAL  FORMULATION 


In  this  Section,  the  theory  and  numerical  algorithms  used  in 
the  structural/hydrodynamic  impact  analysis  are  described.  Since 
the  structural  dynamics  solution  relies  heavily  upon  existing 
software,  the  discussion  of  this  portion  of  the  solution  is  only 
cursory;  for  further  details,  we  encourage  the  reader  to  consult 
the  documentation  for  one  or  more  of  the  structural  analysis 
codes  mentioned  in  Subsection  1.2.  Herein,  emphasis  is  placed 
upon  the  hydrodynamic  solution,  the  finite  elements  used,  and  the 
mesh  rezoning  algorithm. 
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2.1  EQUATIONS  OF  MOTION 

In  what  follows,  a  general  formulation  of  the  soft-body 
impact  problem  is  developed.  Equations  of  motion  are  expressed 


in  terms  of  convected  coordinates  which  coincide  with  initial  or 


reference  positions  of  a  body. 

Consider  two  unconnected  domains  V  and  V,  ,  whose  union  is 

a  b 


denoted  V,  with  boundary  3V  (Figure  2.1).  The  boundary  consists 

of  three  segments:  on  which  external  forces  are  prescribed, 

3v  ,  on  which  displacements  are  specified,  and  aV  =  a V  f)  3V,  , 
u  cab 

over  which  the  two  bodies  are  in  contact.  The  initial  position 

of  a  material  point  in  V  is  described  in  terms  of  Cartesian 

coordinates  a  :m=l,2,3,  and  its  position  at  time  t  is  x.(a  ,t); 

m  i  m 

1=1 , 2 ,3 . 

The  motion  of  both  bodies  must  obey  the  momentum  equation 


(2.1) 


the  traction  boundary  conditions 


on  3  V 


o 


(2.2) 


the  kinematic  boundary  conditions 


u  .  (a  t)  =  u  .  ( t)  on  3V 

i  m ,  l 


u 


(2.3) 


the  contact  interface  conditions 


(2.4) 


» 


c 


(2.5) 


and  the  initial  conditions 


x. (a  ,0) 
i  m 


(2.6) 


(2.7) 
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Figure  2.1.  Contact  Between  Two  Arbitrary  Bodies. 
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t 


In  Equations  ( 2 . 1 ) - ( 2 . 7 ) ,  all  force  and  stress  measures  refer  to 
the  instantaneous  configuration,  and  overbars  signify  prescribed 
quantities. 

We  assume  that  all  kinematic  conditions  (Equations  2.3  and 
2.5)  are  satisfied  identically.  Equations  (2.1),  (2.2)  and  (2.4) 
may  be  expressed  in  the  weak  form 


dA 


3V 


I  [Pjf-Ojfjnj  dA 


(2.8) 


where  fix^^  are  virtual  displacements  which  satisfy  the  kinematic 
boundary  conditions  but  are  otherwise  arbitrary.  The  divergence 
theorem  gives  the  relation 


/  (  o  . .  fix  .  )  .dV  -  /  o..n.6x.  dA  +  /  [  o  .fix.  dA  (2.9) 

y  J1  1  »J  J 1  J  1  9V  J 1  J1  J1 


in  which 


/( o..6x.)  .dV  =  /  o..  .fix.  d  V  +  /  o..6x.  .  dV 


Ji  i  » J 


ji.  j  i 


ji  i  *  j 


Equation  (2.8)  can  therefore  be  rewritten  in  the  more  useful  form 


/  (px.  fix.  + 
V  1  1 


pf.6x.)dV  -  /  t.fix.  dA  -  0 

3V  1  1 

0 


(2.10) 


The  above  weak  form  of  the  momentum  equation  and  force  boundary 
conditions,  together  with  appropriate  constitutive  relationships, 
provides  a  basis  for  the  finite  element  discretization  of  both 
material  regions. 
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2.2  FINITE  ELEMENT  APPROXIMATION 


The  spatial  discretization  of  Equation  (2.10)  is  based  upon 
a  finite  element  approximation  which,  for  a  typical  element  "e", 
has  the  form 


Xi(am  f)  = 
1  ID/  t 


N 

re 

K=1 


*K  (5,n,c) 


1 3 


(t) 


(2.11) 


Here  the  upper  case  subscript  refers  to  individual  nodes  of  the 
element.  The  quantities  x®k(t)  are  nodal  values  of  the  spatial 
coordinates,  which  are  made  continuous  between  adjacent  elements 
using  standard  assembly  procedures  [12].  The  shape  functions 
♦K<5,n,0  are  uniquely  defined  by  the  properties 


W  V‘i> 


in  V 

e 


(2.12) 


♦k  -  0 

Here  is  the  Kronecker  delta, 

equal,  and  equals  zero  otherwise, 
superscript  "e"  without  ambiguity 


outside  V 

e 

which  equals  one  when  I  and  K  are 
In  most  cases,  we  can  omit  the 
,  and  simply  write 


Xi ^am, P  xiK  (2.13) 

Note  that  summation  is  implied  on  the  repeated  index  K,  whose 
range  is  equal  to  the  number  of  nodes  connected  to  an  element. 

The  weak  form  of  the  equations  of  motion  (Equation  2.10) 
becomes,  in  finite  element  form, 


'  [<l>XiJ*jVXiK  +  °ji*K,j4xiK  -  ‘,?iVxiK>  dv  (2-14) 


*  '  .  *i  ♦k5x1K  dA]  '  0 
v  a 


Since  the  nodal  virtual  displacements  6xiK  are  arbitrary  and 
independent  of  one  another,  their  coefficients  must  vanish 
identically.  This  condition  gives  the  semi-discrete  equations  of 
motion  for  the  finite  element  model, 
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E 

£  C  i 

e=l  V 

e 


PXiJ*j*K  +  °ji*k#j 


pf.*K)  dv 


-  I 


3  V 


'KSxiK  dA]  =  0 


(2.15) 


Defining  the  element  mass  matrix 


MJK  -  '  °*j*K  dV 


(2.16) 


and  the  element  force  vector 


F  iK  =  /  (p?i*K-0  ji4,K,jK)dV  +  ;  e  ^i*K  dA  (  2  •  1 7  ) 


Equation  (2.15)  becomes 

\  <HJK  XiJ  -  FtK>  -  0  <2-18) 

e=l 

Finally/  invoking  the  finite  element  assembly  procedure  gives  the 
system  of  equations 


M  x  =  F 
JK  iJ  iK 


(2.19) 


in  which  the  nodal  (upper  case)  subscripts  range  from  one  to  the 
number  of  nodes  in  the  model.  Equations  (2.19)/  subject  to  the 
kinematic  boundary  conditions  (Equations  2.3  and  2.5)  and  the 
initial  conditions  (Equations  2.6  and  2.7)  must  be  integrated  in 
time  to  determine  the  transient  response  of  the  finite  element 
model. 
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2.3  STRUCTURAL  RESPONSE  SOLUTION 


The  structural  responses  of  interest  are  dominated  by  low 
frequency  wave  motion,  and  involve  moderate  amounts  of  material 
deformation.  Compared  with  these  lower  frequency  motions,  the 
highest  frequencies  of  the  computational  mesh  for  the  structure 
are  many  times  greater  due  to  high  material  stiffness.  To 
capture  the  overall  structural  motion  using  time  steps  which  are 
not  controlled  by  high  frequency  behavior  of  the  mesh,  an 
implicit  method  of  solution  is  appropriate.  Implicit  solution 
techniques  make  use  of  the  equations  of  motion  at  time  t+At  to 
solve  for  the  motion  at  t+At,  and  normally  involve  the  solution 
of  simultaneous  equations  and/or  iteration  at  each  time  step. 

As  the  basis  of  the  implicit  integration  we  adopt  the  two 
parameter  Newmark  algorithm  [20],  which  is  obtained  from  the 
following  finite  difference  approximations: 


*  t+At 
x  .  , 

1 J 


t+At-, 
iJ  J 


(2.20) 


t+At 
x  .  , 
iJ 


+ 


+ 


•t  +  At 
x  .  _ 
l  J 


(At)2 


Parameters  a  and  6  are  arbitrary  and  may  be  chosen  to  obtain 
desirable  numerical  properties.  The  values  a=l/4,  6=1/2  give  the 
trapezoidal  integration  rule,  which  is  used  exclusively  in  this 
report . 

In  the  implicit  solution,  Equation  (2.19)  is  applied  at  time 

t+At  to  determine  the  solution  x^+At: 

iJ 


M 


JK 


••  t  +  A  t 
X  .  , 

1 J 


,  t+A  t 
iK 


(2.21) 


Since  the  forces  ,  which  include  the  internal  forces, 

t+At  .  1K 

upon  xiJ  in  a  nonlinear  fashion.  Equation  (2.21)  must  be 


depend 

solved 


12 


iteratively. 


We  define  the  nodal  residual  forces 


t  +  At 
iK 


M 


't  +  At 


JK  i  J 


(2.22) 


which  must  vanish  at  the  true  solution  Zeros  of  Equation 

(2.22)  can  be  generated  by  Newton-Raphson  iteration  based  upon 
the  Taylor  series  approximation 


.  ,-t+Ats 
*^i  K  XMN  ) 


.  ,  t+At,. 

$iK(xMN  ) 


9$ 


iK 


9x 


jL 


( yt+At_  t+Atx 
UjL  XjL  '  * 


0  (2.23) 


The  Jacobian  matrix  which  must  be  formed,  factored  and  solved  to 

carry  out  the  iteration  is  otained  by  substituting  the 

accelerations  x^t^  from  Equations  (2.20)  into  Equation  (2.22), 

1  J 

and  differentiating  with  respect  to  x^.  Notice  that,  since  the 
quantities  x^j  are  discrete  parameters,  we  have  Sx^j/ 9x j 6 ^ j^JL* 
The  result  is: 


9*iK 

9XjL 


9F. 


iK 


1 


3XjL 


a(  At ) 


2  MLK6ij 


(2.24) 


where  6.  .  is  the  Kronecker  delta.  The  leading  term  in 
i  J 

Equation  (2.24)  is  commonly  called  the  tangent  stiffness  matrix. 

Equation  (2.23)  is  applied  repeatedly  at  each  step  of  the 

.  -t  +  At  t  +  At  ,  /  .  . 

dynamic  solution,  until  the  corrections  x^L  ~xjL  and/or  the 

residuals  are  smaller  than  predefined  tolerances.  Once  the 

1 K. 

displacement  solution  is  known  to  the  prescribed  accuracy,  nodal 
velocities  and  accelerations  are  updated  using  Equations  (2.20). 
Reference  [16]  describes  the  implicit  solution  for  the  structural 
response  in  more  detail. 
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2.4  HYDRODYNAMIC  RESPONSE  SOLUTION 


In  contrast  to  the  structural  motion,  the  impact  behavior  of 
the  soft  body  is  hydrodynamic  in  nature  and  involves  very  large 
material  distortions.  Furthermore,  wave  speeds  in  this  portion 
of  the  computational  mesh  are  relatively  small  due  to  low 
material  stiffness.  Here  a  simpler  method  of  solution  is  most 
appropria  te . 

Integration  of  the  equation  of  motion  in  the  hydrodynamic 
mesh  is  performed  using  the  central  difference  approximations 


•t+At/2  _  1  t+At 

iJ  "  St  ,xij 


) 


1 

At 


,*t+At/2 

(xiJ 


t-At/2. 

iJ  ‘ 


(2.25) 


The  solution  in  this  case  is  explicit,  since  the  state  at  time 
t+At  is  determined  directly  from  that  at  time  t.  Equations 
(2.19)  and  (2.25)  are  used  directly,  once  per  time  step,  to 
advance  the  acceleration,  velocity  and  displacement  solutions. 

The  explicit  nature  of  Equation  (2.25)  makes  the  integration 
conditionally  stable  (that  is,  stable  only  for  a  limited  range  of 
step  sizes).  The  stability  limit  is  well-known  [21],  and  is  given 
by 

Atcr  =  -  lAT?  -  e]  (2.26) 

ma  x 

where  u>  is  the  highest  frequency  of  the  mesh,  and  e  is  the 

ma  X 

fraction  of  critical  damping  in  the  highest  frequency.  Since  it 
is  inconvenient  to  estimate  the  maximum  frequency  of  the  mesh 
directly,  we  note  that  the  maximum  frequency  for  the  complete 
mesh  is  bounded  above  by  the  maximum  element  frequency.  It  is 
therefore  possible  to  relate  the  highest  frequency  to  the  nodal 
separation  At  and  the  material  sound  speed  c, 
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(2.27) 


Equation  (2.27)  can  be  applied  element-by-element  at  each  time 
step  to  control  the  integration  time  step.  The  acoustic  speed  c 
for  an  inviscid  flow  is  simply 

c  =  |^)  (2.28) 
'  s 

where  the  symbol  s  represents  the  entropy.  This  quantity  can  be 
obtained  for  each  element  as  a  by-product  of  the  equation  of  state 
calculations. 
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2.5  HYDRODYNAMIC  FINITE  ELEMENTS 


The  details  of  the  specific  finite  element  approximation 
adopted  for  the  hydrodynamic  solution  are  described  in  this 
section.  Basic  aspects  of  the  element  formulation  follow  the  work 
of  Flanagan  and  Belytschko  [21]. 

2.5.1  Element  Geometry  and  Shape  Functions 

The  unit  of  approximation  in  the  hydrodynamic  mesh  is  the 
eight-node,  linear  displacement  hexahedron  shown  in  Figure  2.2.  A 
unit  cube  in  the  parametric  coordinates  (€,n,c),  also  denoted  by 
is  mapped  into  physical  coordinates  x^  by  means  of  the  shape 
functions  If  the  center  of  the  cube  is  taken  to  be 

the  origin,  -l/2_<C^l/2 ,  the  shape  functions  are 

♦  jU/TWC)  =  (§+2?I5)  (■j+2nIn)  (^+2CjC)  (2.29) 

Flanagan  and  Belytschko  [21]  have  shown  that  the  <|> x  also  may 
be  expressed  in  terms  of  an  orthogonal  set  of  base  vectors 


*1  8ZI  +  45A1I  +  4nA2I  +  4U3I 


(2.30) 


+  +  7CCr2l  +  ?^r3I  +  CnCr4I 


The  constants  /  and  raI  are  summarized  in  Table  2.1,  and 

are  hereafter  referred  to  as  the  constant,  linear,  and  hourglass 
base  vectors,  respectively.  This  form  of  the  shape  functions  is 
useful  since  it  isolates  the  contribution  of  hourglass  modes  of 
deformation  (which  are  neglected  by  one-point  integration) 
explicitly  on  the  unit  cube. 

2.5.2  Mean  Stress  Approximation 

The  eight-node  hexahedron,  one  of  the  simplest  of  three 
dimensional  elements,  must  be  used  with  care  if  good  numerical 
behavior  is  to  be  obtained.  An  exact  numerical  integration  of 
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4 


Figure  2.2.  Eight-Node  Hexahedral  Finite  Element. 
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TABLE  2.1 

BASIS  VECTORS  FOR  EIGHT-NODE  HEXAHEDRON 
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the  eight  node  element  may  lead  to  "locking"  of  the  mesh,  while  a 
simple  reduced  (one  point)  integration  fails  to  integrate  the 
volume  of  irregular  elements  exactly  and  may  be  non-convergent. 

In  the  present  work,  we  adopt  the  "mean  stress  approximation" 
described  in  Reference  [21]  to  obtain  a  convergent  element  with 
desirable  numerical  properties. 

The  mean  stress  approximation  is  obtained  by  using  mean 
values  of  the  stress  and  kinematic  variables,  while  performing 
the  volume  integration  exactly  for  arbitrarily-shaped  elements. 
The  resulting  element  is  therefore  capable  of  representing  states 
of  uniform  stress  and  deformation  exactly,  a  property  which  is 
necessary  for  convergence  [12].  The  sole  disadvantage  is  that 
hourglass  deformation  modes  of  the  element  may  be  unrestrained, 
and  recourse  must  be  made  to  artificial  damping  mechanisms  to 
avoid  contamination  of  the  numerical  solution. 

2.5.3  Element  Kinematics 

The  spatial  coordinates  of  points  within  a  single  element 
are  approximated  as  indicated  in  Equation  (2.11), 


x 


i 


4*KxiK 


(2.31) 


where  the  shape  functions  for  the  eight  node  hexahedron  are 
defined  by  Equation  (2.30).  The  range  of  the  nodal  (upper  case) 
subscript  is  8,  the  number  of  nodes  per  element.  The  velocity 
field  is  expressed  in  a  similar  form 


vi  =  *KviK 


(2.32) 


For  the  spatial  velocity  gradient,  we  obtain 


v .  .  =  <fc  .  v . 

1/3  ?K,jviK 


(2.33) 


since  only  <t>  varies  with  spatial  position.  The  parameters  v. 


iK 
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are  nodal  quantities  which  vary  only  with  time/  and  hence 
ViK,j  * 

With  the  mean  stress  approximation,  it  is  necessary  to  form 
mean  values  of  the  velocity  gradient  and  stress  for  an  element. 

For  this  purpose,  it  is  convenient  to  define  a  matrix  containing 
the  integrals  of  the  shape  function  derivatives, 

BiJ  =  v/* J,i  dV  (2*34) 


For  lack  of  a  better  name,  matrix  B  will  be  referred  to  as  the 
"geometric  matrix"  for  an  element,  since  it  is  purely  a  function 
of  the  element  nodal  coordinates.  Consider  the  mean  velocity 
gradient  for  an  element, 


v . 
i/l 


(2.35) 


From  Equation  (2.34), 


v .  . 

i/l 


can  be  written  more  concisely  as 


v .  . 

i/l 


1 

V 


iK 


(2.36) 


The  geometric  matrix  B  also  provides  a  convenient  method  for 
evaluating  the  element  volume,  since  [21] 


x 


iKB  jK 


(2.37) 


where  is  the  Kronecker  delta.  It  is 
that  both  the  constant  and  the  hourglass 
orthogonal  to  B, 


straightforward  to  verify 
base  vectors  are 


BjIzI  -  0 

7  j  — 1/2/3 

(2.38) 

Bjirax=  0 

7  a=l / 2 , 3 / 4 

(2.39) 

20 


2.5.4  Constitutive  Model 


The  constitutive  relations  used  for  the  hydrodynamic  mesh 
are  relatively  simple,  since  the  larger-scale  effects  of  impact 
(such  as  pressure  level  and  momentum  transfer)  are  of  primary 
interest.  The  material  is  defined  in  terms  of  bulk  and  shear 
moduli  ( K , G ) ,  density  p,  and  a  material  strength  a  .  The  rate  of 
deformation  and  the  co-rotationa 1  (Jaumann)  rate  of  the  Cauchy 
stress  are  assumed  to  be  related  by 


a.  .  =  (K-|g)  d,  .  6  .  .  +  2G  d.  . 
13  3  kk  13  13 


(2.40) 


where  the  Jaumann  stress  rate  is  defined  by 


a.  .  =  5.  .  +  a.T, w.  .  +  a.,  w,  . 
ij  13  lK  k}  jk  ki 


(2.41) 


In  Equation  (2.41),  a.,  is  the  stress  tensor,  o.. 
time  derivative,  and  is  the  spin  (vorticity). 


its  material 


v . 

1/ 


j 


(2.42) 


Material  strength  is  accounted  for  by  computing  an  effective 
deviatoric  stress. 


o 

e 


o!  .o'.  . 
1 3  13 


(2.43) 


which  is  limited  in  magnitude  to  the  material  strength  a  .  when 
the  strength  is  exceeded,  the  deviatoric  stresses  are  scaled 
according  to 


(2.44) 


It  should  be  noted  that  the  primary  purpose  of  the  viscous 
terms  in  the  constitutive  model  is  to  stabilize  pure  shear  modes 
which  are  otherwise  unconstrained.  For  this  purpose,  a  small 
positive  shear  coefficient  is  normally  used. 
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2.5.5  Element  Forces 

The  element  forces  (Equation  2.19)  include  three  separate 
contributions: 

-  prescribed  external  forces: 

-  internal  forces;  and 

-  hourglass-stabilizing  forces. 

The  first  two  of  these  are  physically  meaningful.  The  hourglass 
resisting  forces  exist  solely  to  stabilize  the  computational  mesh 
and  must  not  interfere  with  the  representation  of  the  problem 
physics . 

Physical  forces  for  an  element  are  obtained  from  Equation 
(2.17).  For  simplicity,  we  will  consider  only  body  forces  and 
the  element  internal  forces, 


F 


ik 


/ 

V 


(  Pf  i  * 


K 


a  .  . 
Dl 


dV 


(2.45) 


Consistent  with  the  mean  stress  approximation, 
replaced  by  their  mean  values  over  an  element; 
nodal  forces  are  then 


p ,  f . ,  and  a  .  . 
i  li 

the  resulting 


a  re 


iK 


=  8  PVVK 


o  .  .  B  . 
11  IK 


(2.46) 


The  mean  stress  o__  is  obtained  directly  from 
model  and  the  mean  velocity  gradient. 


the  constitutive 


Hourglass-resisting  forces  for  the  hexahedron  are  computed 
using  the  method  described  by  Flanagan  and  Belytschko  [21]. 
Figure  2.3  shows  the  four  independent  displacement  modes  of  an 
element  corresponding  to  a  single  component  of  displacement;  the 
purpose  of  the  anti-hourglassing  forces  is  to  resist  the 
development  of  the  four  hourglass  displacement  modes,  which  are 
strain-free  under  the  mean  stress  approximation. 
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Figure  2.3.  Hourglass  Displacement  Patterns  for  Eight-Node 
Hexahedron . 
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as 


Define  the  hourqlass  shape  vector  y 

1  „  aI 


y  t  =  r„T  -  T7  B  •  -r x  . , r  , 
“I  “I  V  ll  iJ  aj 


(2.47) 


The  first  subscript  refers  to  a  particular  hourglass  displacement 
mode  (Figure  2.3) ,  and  ranges  from  one  to  four;  the  upper  case 
subscript  refers  to  specific  nodes  of  an  element.  The  hourglass 
shape  vector  is  orthogonal  to  the  linearly  varying  portion  of  the 
element  velocity  field  for  an  element  of  arbitrary  shape.  That 
is,  the  products 


(2.48) 


called  the  hourglass  modal  velocities/  provide  an  indication  of 
the  amount  of  hourglassing  which  is  present  in  the  element.  The 
factor  of  1//8  in  Equation  (2.48)  is  present  to  normalize  the 
magnitudes  of  the  modal  velocities.  Anti-hourglassing  forces  are 
proportional  to  the  hourglass  modal  velocities/ 


(2.49) 


and  vanish  when  the  element  velocity  field  is  purely  linear.  A 
damping  coefficient  e  is  used  to  control  the  magnitude  of  the 
restoring  forces.  When  the  anti-hourglassing  forces  are  defined 
as  indicated  in  Equation  (2.49)/  e  corresponds  to  the  fraction  of 
critical  damping  in  the  highest-f requency  mode  of  the  element/ 
and  critical  time  steps  may  be  estimated  using  Equation  (2.26). 

It  is  important  to  realize  that  only  pure  hourglassing 
deformations  are  suppressed  by  the  above  scheme.  When  global 
deformation  modes  develop  in  which  individual  finite  elements 
experience  other  than  purely  linear  velocity  fields,  the  anti- 
hourglassing  forces  combine  to  produce  zero  resultant  forces, 
since  the  hourglass  modal  velocities  in  adjacent  elements  are 
identical  in  sign.  In  contrast  to  global  deformation  modes,  true 
hourglassing  patterns  contain  modal  velocities  which  alternate  in 
sign  between  adjacent  elements,  and  the  resistive  forces  combine 
to  suppress  such  motion  (Figure  2.4). 
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Figure  2.4.  Behavior  of  Hourglass-Resisting  Forces  in 
Adjacent  Finite  Elements. 
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2.6  MESH  REZONING 


In  the  hydrodynamic  portion  of  the  mesh,  large  deformations 
may  occur  as  the  solution  proceeds.  Since  the  nodal  positions 
are  updated  at  each  time  step,  very  large  deformations  may  result 
in  excessive  distortion  of  the  mesh,  which  are  manifested  in  the 
following  ways: 

-  disparities  in  element  size  due  to  volume  change; 

-  singularities  due  to  element  volumes  tending  to  zero; 

-  undesirable  element  shape  due  to  large  distortions;  and 

-  reductions  in  the  permissible  time  step  due  to  extreme 

distortion  and/or  volume  change. 

These  unfortunate  situations  can  be  resolved  by  reconstruction  of 
the  computational  mesh  (rezoning) ,  performed  as  necessary  during 
the  solution. 

Rezoning  techniques  are  commonly  used  in  large  deformation 
problems  of  solid  mechanics,  such  as  ballistic  impact  and  metal 
forming.  However,  their  use  normally  requires  a  great  deal  of 
interaction  with  the  analyst  and  is  therefore  time-consuming  and 
expensive.  The  rezoning  method  developed  in  the  present  work  is 
automatic  and  non-interactive,  and  therefore  can  be  performed 
without  halting  the  solution  temporarily.  In  return  for  this 
convenience,  our  technique  allows  somewhat  less  control  over  the 
appearance  of  the  rezoned  computational  mesh. 

The  rezoning  process  involves  redefining  nodal  coordinates 
and  element  connections  for  the  model,  and  relocating  solution 
variables  (velocities,  pressures,  densities,  stresses)  at  the 
newly-defined  nodal  positions.  In  general,  the  size  of  the 
rezoned  model  (number  of  nodes  and  elements)  may  be  different 
from  that  of  the  original  model.  Subsequent  sections  describe 
the  conditions  under  which  mesh  rezoning  is  performed,  and  the 
main  features  of  the  rezoning  algorithms  used  in  the  present 
work . 

Two  types  of  rezoning  procedures  have  been  developed  and  used 
to  date;  we  refer  to  these  as  rezoning  algorithms  I  and  II  in  the 
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remainder  of  this  section.  The  logic  for  invoking  the  rezoning 
process  is  identical  for  both  algorithms ,  and  is  discussed  in  the 
next  subsection. 

2.6.1  Initiation  of  Rezoning  Procedure 

In  most  instances,  severe  distortion  of  the  mesh  can  be 
detected  by  monitoring  element  volumes  and  computed  critical  time 
step  sizes.  If  uniform  compression  or  expansion  occurs,  critical 
time  step  sizes  may  change  slowly  while  the  elemental  volumes 
change  more  rapidly.  Severe  shape  distortion  is  usually 
accompanied  by  a  rapid  decrease  in  the  allowable  time  increment. 

The  hydrodynamic  solution  module  continuously  monitors  both 
the  minimum  and  maximum  element  volumes  for  the  entire  mesh,  as 
well  as  the  allowable  time  step.  Each  time  the  hydrodynamic 
solution  is  restarted,  initial  values  of  these  parameters  are 
computed  and  stored.  At  each  time  step,  current  values  are 
recovered  from  the  internal  force  calculation,  and  changes  in 
each  quantity  are  compared  with  preset  tolerances  (user-specified 
or  defaulted) .  When  any  one  of  the  monitored  values  (minimum 
volume,  maximum  volume,  critical  time  step),  as  compared  with  the 
initial  values,  exceeds  its  tolerance,  rezoning  is  automatically 
invoked.  Following  the  rezone,  several  key  solution  quantities 
(e.g.,  nodal  masses)  and  the  initial  volume  and  time  step  values 
are  reset,  and  the  solution  restarts  from  the  point  of 
interruption. 

2.6.2  Rezoning  Algorithm  I 

The  first  rezoning  procedure  consists  of  two  phases:  (a) 
the  definition  of  new  nodal  points  and  the  corresponding  solution 
variables;  and  (b)  the  definition  of  finite  elements  which  span 
the  material  volume.  The  algorithm  is  designed  to  produce  a 
regular  mesh  within  the  interior  of  the  region,  with  irregular 
elements  used  as  necessary  near  the  boundaries  and  free  surfaces 
of  the  model. 

The  relocation  of  nodal  points  within  the  mesh  is  performed 
using  a  reference  grid  pattern  specifying  fixed  mesh  stations 
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along  each  of  the  Cartesian  coordinate  axes.  Figure  2.5  shows 
such  a  reference  grid  pattern  in  two  dimensions.  Wherever 
material  points  coincide  with  intersections  on  this  reference 
pattern,  nodes  will  be  positioned;  on  the  boundaries  of  the 
region,  nodes  are  located  on  reference  lines  as  much  as  possible. 

Figure  2.6  shows  the  relocation  procedure  in  a  conceptual 
form.  All  nodes  within  the  interior  of  the  body  lie  at  reference 
grid  intersections,  and  additional  nodes  (on  reference  lines  if 
possible)  describe  the  boundary  of  the  region  as  required.  As 
the  positioning  of  new  nodes  takes  place,  linear  interpolations 
within  the  existing  finite  elements  are  performed  to  define  nodal 
values  of  the  solution  variables  at  the  new  node  locations. 

A  key  ingredient  of  the  numerical  algorithm  used  in  rezoning 
is  the  establishment  of  a  one  to  one  correspondence  between  node 
points  and  mesh  reference  positions.  When  a  conflict  exists,  it 
is  resolved  in  favor  of  the  (new)  node  which  lies  closest  to  the 
boundary  of  the  material.  This  convention  facilitates  searching 
operations  during  nodal  placement,  the  merging  of  boundary  nodes, 
and  the  creation  of  finite  elements  for  the  rezoned  mesh.  In 
return  for  this  convenience,  the  resolution  of  the  procedure  is 
limited  to  a  linear  interpolation  between  points  which  lie  on 
adjacent  reference  lines. 

This  first  phase  of  the  rezoning  process  is  accomplished  by 
examining  each  element  once,  determining  its  intersections  with 
the  reference  mesh  pattern  and  then  tentatively  positioning  new 
nodal  points.  Although  new  nodes  may  be  moved  several  times  as 
the  rezone  continues,  this  operation  is  quick  since  no  searching 
is  required.  Once  all  elements  have  been  searched,  the  boundary 
nodes  are  examined  for  excessively  small  spacing  and  merged  when 
possible;  this  prevents  critical  time  step  sizes  for  the  solution 
phase  from  becoming  unreasonably  small  due  to  fine  mesh  spacing 
where  none  is  called  for. 

Creation  of  new  elements  is  straightforward  once  the  nodal 
points  have  been  redefined  and  merged.  Each  set  of  eight 
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Figure  2.5.  Typical  Reference  Grid  for  Mesh  Rezoning. 
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Figure  2.6.  Two-Dimensional  Example  of  Rezoning  Procedure. 
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adjacent  reference  positions  potentially  defines  one  element;  to 
determine  if  an  element  is  actually  present,  the  nodes  which 
correspond  to  each  such  group  of  points  is  examined.  When  at 
least  four  nodes  are  present,  and  the  resulting  element  volume  is 
positive,  an  element  is  defined. 

2.6.3  Rezoning  Algorithm  II 

A  second  rezoning  technique  has  been  developed  based  upon 
experiences  gained  in  using  Algorithm  I.  The  function  of  this 
alternative  technique  is  the  same:  given  a  distorted  mesh  of 
eight-node  elements,  generate  a  revised  mesh  in  which  the 
distortion  is  reduced  or  eliminated. 

Algorithm  II  is  based  upon  a  different  philosophy  from 
Algorithm  I,  in  that: 

o  the  number  of  nodes  and  elements  does  not  change; 

o  mass  and  momentum  are  conserved  to  machine  accuracy;  and 

o  the  repositioning  is  less  dramatic  than  with  Algorithm  I. 

This  second  procedure  is  performed  by  examining  each  node  of  the 
existing  mesh,  and  determining  whether  the  node  lies  within  the 
interior  of  the  body,  on  a  smooth  surface,  on  a  distinct  "edge"  of 
the  region,  or  at  a  current  "corner."  The  distinction  between 
these  cases  is  made  by  examining  the  element  faces  connected  to 
each  node. 

For  each  point  in  the  model,  all  element  faces  connected  to 
the  point  are  collected,  and  a  geometric  analysis  much  like  a 
hidden  line/surface  removal  procedure  is  performed.  For  an 
interior  node,  all  connected  surfaces  are  eliminated  (thus 
identifying  the  node  as  interior) .  In  this  case,  neighboring 
points  are  used  to  reposition  the  node  in  question,  and  solution 
variables  are  mapped  to  the  new  geometry  by  conditions  of  mass  and 
momentum  conservation,  and  by  volume  weighting  of  tensor-valued 
quantities  (such  as  stress) .  A  series  of  similar  tests  have  been 
designed  to  isolate  surface,  edge,  and  corner  nodes.  The  sole 
difference  in  treating  these  additional  cases  lies  in  the 
selection  of  surrounding  nodes  to  be  used  in  repositioning.  For 
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surface  nodes,  only  those  nodes  on  the  same  surface  are  used,  in 
order  to  maintain  the  surface  shape;  "edge”  nodes  are  moved  only 
along  the  curve  defining  the  edge  in  question.  Corner  points  are 
never  relocated,  since  they  represent  locations  at  which  the  model 
geometry  may  be  sharply  discontinuous. 

Our  experience  has  been  that  the  second  rezoning  algorithm  is 
superior,  since  mass,  momentum,  and  geometric  shape  are  preserved. 
Symmetry  conditions,  when  they  exist,  tend  to  survive  many 
repeated  rezones  with  Algorithm  II  as  well.  Finally,  the  second 
algorithm  requires  somewhat  less  memory  and  is  faster  in  execution 
than  Algorithm  I. 

The  sole  disadvantage  to  the  alternative  rezoning  procedure 
is  that  nodal  repositioning  is  moderate  compared  with  the  first 
algorithm,  so  that  rezoning  must  be  performed  more  frequently. 

The  most  satisfactory  performance  has  been  obtained  with 
relatively  mild  tolerances  (about  50%)  for  both  maximum  element 
volume  change  and  time  step  reduction. 

The  soft-body  impact  code  has  been  written  to  use  either  of 
the  rezoning  algorithms  described  in  this  section.  However, 
Algorithm  II  is  now  being  used  exclusively  due  to  its  superior 
performance . 
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2.7  STRUCTURE  MODEL  /  HYDRODYNAMIC  MODEL  INTERFACE 

The  interface  between  structural  and  hydrodynamic  models 
uses  the  contact  analysis  technique  described  in  References  [16] 
and  [18].  "Interface"  or  "contact"  elements  are  defined  on  the 
surface  of  the  structure  model:  these  elements  are  strictly 
geometric,  and  have  no  material  characteristics.  The  interface 
elements  move  with  the  structure  model,  and  cannot  be  penetrated 
by  nodes  of  the  hydrodynamic  model.  Thus,  the  interface  elements 
impose  time-dependent  displacement  boundary  conditions  on  the 
hydrodynamic  model.  Forces  at  hydrodynamic  nodes  in  contact  with 
the  interface  elements  define  the  impact  forces,  which  are 
transferred  to  the  structural  model  at  each  major  time  step.  At 
minor  time  steps  (i.e.,  those  in  the  soft-body  mesh),  only 
kinematic  calculations  (related  to  the  displacement  boundary 
conditions)  are  performed. 

The  interface  elements  used  in  the  present  work  are  defined 
by  four  corner  points,  as  shown  in  Figure  2.7.  Ordering  of  the 
nodal  points  determines  the  outward  normal  direction  for  an 
element.  The  restriction  to  four  points  rather  than  the  variable 
number  of  nodes  used  in  the  original  algorithm  is  for  reasons  of 
computational  economy,  since  contact  calculations  must  be 
performed  at  each  time  step  for  the  hydrodynamic  mesh. 

The  sole  remaining  differences  between  the  interface  element 
calculations  in  the  present  work  and  the  original  algorithm  are 
related  to  the  ordering  of  the  computations.  Since  the  interface 
element  positions  change  only  at  major  (structure)  time  steps, 
while  contact  searches  must  be  done  at  minor  (hydrodynamic)  time 
steps,  geometric  parameters  for  the  interface  are  precomputed  at 
each  major  step  and  stored  for  repeated  use.  These  parameters 
include  element  local  coordinate  transformations,  coordinates,  and 
limiting  coordinate  values. 
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Figure  2.7.  Structural-Hydrodynamic  Interface  Element. 
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SECTION  3 

DEMONSTRATION  PROBLEMS 


Testing  of  the  fluid-structure  impact  analysis  procedure  has 
included  the  following  types  of  problem  solutions: 

o  elastic  vibrations  (using  the  moving  boundary  option); 
o  low-speed  normal  impacts  on  rigid  surfaces; 
o  classical  fluid  flow  cases,  such  as  Couette  flow  (using  a 
Newtonian  fluid  material  model); 
o  oblique  impacts  on  rigid  surfaces;  and 
o  coupled  solutions  for  normal  fluid-structure  impact. 

In  most  cases,  these  solutions  have  been  used  to  verify  the 
correctness  of  the  soft-body  impact  analysis  module.  However, 
some  cases  (mostly  oblique  Impacts)  have  been  used  only  to  study 
the  behavior  of  the  ant i-hourglassing  formulation  and  to  verify 
the  correct  operation  of  the  rezoning  module  for  problems 
involving  fully  three-dimensional  response.  A  fully  coupled 
impact  test  case  is  described  in  detail  below. 

The  coupled  problem  considered  is  the  normal  impact  of  a 
cylinder  of  porous  gelatin,  used  in  the  laboratory  simulation  of 
bird  impacts,  on  a  square  flat  plate.  Although  experimental  data 
corresponding  precisely  to  the  problem  considered  do  not  exist,  a 
similar  case  of  impact  on  a  rigid  target  is  reported  in  Reference 
[233.  The  objectives  in  studying  a  similar  case  with  a  flexible 
target  are  to  verify  that  the  computed  force  levels  are  similar, 
and  to  observe  the  performance  of  the  contact/impact  procedure  in 
a  situation  where  kinematic  constraints  and  interface  forces  can 
be  observed  and  interpreted  simply. 

Physical  data  used  in  the  simulation  are  as  follows.  The 
plate  is  12  inches  square  and  0.1  inch  thick;  one  quadrant  of  the 
plate  is  modeled  due  to  symmetry.  The  plate  material  is 
aluminum,  with  a  modulus  of  10  Mpsi,  Poisson’s  ratio  0.3.  and 
density  0.098  pounds  per  cubic  inch.  All  of  the  outer  boundaries 
of  the  plate  are  fully  clamped.  For  simplicity,  and  to  preclude 
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the  development  of  very  large  deformations,  the  plate  material  is 
assumed  to  behave  elastically;  however,  large  displacement 
effects  are  included  in  the  computation.  The  impacting  body  is  a 
cylinder  1  inch  in  radius  and  4  inches  long  (L/D=2).  The 
cylinder  material  is  gelatin  with  10*  porosity,  with  density 
0.0344  pounds  per  cubic  inch.  No  mechanical  properties  are 
available  from  Reference  [23];  for  this  calculation,  we  assume  a 
shear  modulus  of  100.  psi,  linear  and  quadratic  bulk  moduli  of 
1000.  psi,  and  a  material  strength  of  10.  psi. 

The  impact  velocity  for  the  event  has  been  chosen  to 
correspond  to  one  of  the  rigid  target  cases  reported  in  Reference 
[23],  at  4724.4  inches  per  second  (120  meters/second).  The 
approximate  "squash-up"  time  (the  time  for  the  cylinder  to  travel 
a  distance  equal  to  its  initial  length)  is  therefore  0.0008466 
seconds . 

For  the  finite  element  solution  of  this  problem,  25  thin 
shell  elements  are  used  to  represent  the  plate,  and  48  eight-node 
hexahedra  model  the  cylinder.  Both  of  these  models  are  coarse, 
since  the  objective  is  not  to  perform  a  detailed  stress  analysis 
of  the  event;  the  most  important  effects  which  the  solution  must 
capture  are  the  transfer  of  momentum  between  cylinder  and  plate, 
and  the  tendency  of  the  impacting  body  to  spread  over  the  surface 
of  the  target.  In  the  structure  mesh,  a  constant  time  step  of 
0.01  milliseconds  (about  1.2*  of  the  squash-up  time)  is  used;  for 
the  soft  body  mesh,  the  time  step  is  adjusted  automatically  and 
continuously  within  the  program.  It  is  important  to  notice  that 
the  resolution  with  which  impact  forces  can  be  obtained  from  the 
solution  is  limited  to  the  structure  mesh  time  step,  since  the 
interface  force  calculations  are  performed  only  at  the  beginning 
of  each  time  step. 

Figure  3.1  gives  an  overview  of  the  impact  calculation, 
showing  the  boundaries  of  both  bodies  at  times  of  0.1  ms,  0.3  ms, 
0.5  ms,  and  0.7  ms.  The  configuration  at  0.7  ms  represents  the 
final  point  of  the  solution,  at  which  time  the  cylinder  just 
begins  to  rebound  from  the  surface  of  the  plate.  By  this  time, 
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Geometry  of  Plate  and  Cylinder  at 
0.3  ms,  0 . 5  ms ,  and  0.7  ms . 


the  plate  has  attained  its  maximum  deflection  and  is  vibrating 
upward.  Obviously  this  same  rebounding  does  not  occur  for  the 
rigid  plate  tests  reported  in  Reference  [23]. 

Figure  3*2  shows  the  deformed  geometry  of  the  cylinder  at 
intervals  of  0.05  ms.  The  body  is  about  70 %  consumed  during  the 
solution,  as  shown  in  the  side  view  of  Figure  3.3  (just  at  the 
point  of  rebound).  This  is  at  least  qualitatively  correct:  the 
deflection  of  the  plate  reaches  about  0.75  inches  at  0.65  ms, 
which  means  that  the  rear  surface  of  the  cylinder,  traveling  at 
its  initial  velocity,  would  reach  the  plate  at  approximately  1.0 
ms.  The  soft-body  mesh  has  been  rezoned  five  times  during  the 
interval  shown  in  the  Figure.  The  usual  cause  for  rezoning  is 
exceedance  of  the  time-step-change  tolerance,  which  is  set  to 
0.75;  in  two  instances,  the  element  volume  change  criterion  has 
been  activated  due  to  extreme  compression  of  elements  adjacent  to 
the  impact  surface. 

Time  histories  of  the  interface  forces,  vertical  momentum  of 
the  cylinder,  and  plate  central  deflection  are  given  in  Figures 
3.4  through  3.6,  respectively.  The  theoretical  impact  force 
during  steady  flow,  for  the  quarter  of  the  cylinder  modeled,  is 
[233: 


F  =  1  pAuQ2  =  1558.4  lb. 

This  force  level  is  indicated  by  a  horizontal  line  in  Figure  3.4. 
During  the  initial  stages  of  the  impact,  the  total  vertical  force 
oscillates  about  the  theoretical  value,  and  becomes  larger  at 
later  times.  The  following  important  points  should  be  noted: 

o  the  total  forces  plotted  in  Figure  3.4  represent  values  at 
the  beginning  of  structure  time  steps  only,  and  therefore 
are  limited  in  resolution  (note  that  the  momentum  curve 
for  the  cylinder,  Figure  3.5,  is  much  smoother,  since 
this  output  is  generated  several  times  per  increment); 
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Figure  3.2.  Soft-Body  Mesh  at  Time  Intervals  of 
0.  05  Milliseconds. 
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Figure  3.3.  Elevation  View  of  Cylindrical  Body  at 
0.7  Milliseconds. 
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Figure  3.4.  Impact  Force  versus  Time. 
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Figure  3.5.  Vertical  Momentum  of  Soft  Body  Mesh  versus  Time 
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TIME  (milliseconds) 


o  the  shell  elements  used  for  the  plate  are  bilinear  and 
give  a  rather  coarse  displacement  shape;  as  individual 
elements  change  their  orientation  abruptly  from  step  to 
step,  unrealistic  oscillations  occur  in  the  vertical 
forces  which  would  not  be  expected  with  higher-order 
structural  finite  elements; 

o  the  increase  in  vertical  forces  late  in  the  event  is 
directly  related  to  the  displacement  behavior  of  the 
plate;  note,  for  example,  that  the  peaks  in  the  force 
curve  correspond  to  times  at  which  the  plate  velocity  has 
changed  sign,  and  the  center  of  the  plate  moves  upward. 

As  Figure  3.6  shows,  the  plate  central  deflection  reaches  a 
maximum  near  0.65  ms,  and  exhibits  some  higher-frequency 
oscillations  prior  to  that  time.  These  oscillations  are  rather 
mild  during  the  early  stages  of  the  impact,  but  become  more 
pronounced  and  shorter  in  duration  as  the  deflections  increase, 
presumably  due  to  membrane  stiffening  effects. 

The  results  of  the  normal  impact  case  above  are  encouraging, 
since  the  calculations  are  well-behaved  and  yield  reasonably 
accurate  information  about  the  impact  interface  forces.  Further 
testing  needs  to  be  performed  to  address  the  following  questions: 

(1)  Does  a  finer  (or  higher-order)  structural  model  lead  to 
improved  impact  force  predictions? 

(2)  Does  the  solution  predict  accurately  the  distribution  of 
loading  over  the  target?  Useful  experimental  results  exist  for 
aircraft  transparencies,  with  loading  areas  recorded  using  high 
speed  photography,  which  can  be  used  to  verify  this  aspect  of  the 
analysis  technique. 

(3)  How  accurately  must  the  bird  (or  bird  substitute)  material  be 
characterized  in  order  to  obtain  reasonable  predictions  of  impact 
pressures,  spreading,  and  load-response  coupling? 

We  feel  that  the  most  serious  limitations  which  exist  at 
present  have  not  to  do  with  the  analytical  model,  but  with  the 
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lack  of  reliable  data  for  materials  characterization.  While 
mechanical  property  values  may  not  be  critical  to  the  calculation 
of  steady  flow  phase  impact  forces,  inaccurate  material  modeling 
will  affect  the  spatial  distribution  of  loading  on  the  target, 
particularly  for  more  complex  geometries.  In  many  transparency 
impact  situations,  the  interaction  between  structural  response 
and  the  spatial  distribution  of  the  loading  is  felt  to  be  a 
critical  factor  in  determining  the  transparency  dynamic  response. 
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SECTION  4 

COMPUTER  CODE  DESCRIPTION 


This  Section  contains  a  brief  description  of  the  computer 
implementation  of  the  soft-body  impact  analysis.  Information 
contained  in  the  section  will  be  of  interest  to  programmers  or 
other  persons  involved  in  further  development/  conversions  of  the 
computer  code  to  other  machines/  or  interfacing  of  the  analysis 
routines  with  other  structural  dynamics  programs.  Information 
about  the  routine  usage  of  the  computer  program  can  be  found  in 
Section  5. 
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4.1  OVERVIEW 


The  computer  code  which  implements  the  soft-body  impact 
analysis  consists  of  a  set  of  input/  analysis  and  output 
subroutines/  which  can  be  executed  either  as  a  stand-alone 
program  (with  a  suitable  main  program  inserted)/  or  coupled  with 
existing  programs  for  structural  dynamic  analysis.  For 
convenience/  we  refer  to  this  collection  of  software  as  the  SBI 
(Soft— Body  Impact)  code.  The  existing  finite  element  program 
MAGNA  [16]  has  been  used  as  the  structural  analysis  module/  it  is 
simply  referred  to  by  name  in  the  discussion  to  follow. 

The  SBI  code  is  written  in  ANSI  FORTRAN  77,  with  the 
exception  of  one  subprogram  (CPUSEC)  which  must  be  provided  for 
each  target  machine.  The  program  has  been  compiled  and  executed 
on  CDC ,  CRAY  and  VAX  computers  without  incident.  The  machine 
dependencies  in  SBI  are  limited  to  three  categories: 

(1)  the  function  subprogram  CPUSEC,  which  returns  the  CPU 
time  from  the  start  of  the  job,  is  machine-specifier 

(2)  IMPLICIT  DOUBLE  PRECISION  statements  are  included  in 
each  subprogram  which  uses  variables  of  real  data 
typer  these  IMPLICIT  statements  contain  the  leading 
characters  "C-DBL " ,  which  can  be  removed  to  activate 
the  statements  on  machines  with  short  word  lengthsr 
and 

(3)  file  OPEN  statements  used  in  a  driver  program  in  the 
stand-alone  mode  may  require  different  file  name 
parameters  for  the  printed  output  file  (e.g.,  OUTPUT 
on  CDC,  $OUT  on  CRAY,  TT:  or  none  on  VAX-11). 

The  remaining  subsections  describe  the  major  components  of 
the  SBI  code,  and  the  routines  which  are  necessary  to  interface 
the  hydrodynamics  analysis  with  a  structural  analysis  code. 
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4.2  MODULE  AND  COMMON  BLOCK  SUMMARIES 


The  hydrodynamics  (SBI)  code  consists  of  five  primary 
program  modules:  input/  analysis/  rezoning/  output  and  restart. 
These  modules  can  be  invoked  in  virtually  any  order/  except  that 
the  input  module  must  be  called  once  at  the  beginning  of  a  new 
problem.  Entry  points  into  the  first  three  modules  are  unique: 
that  is/  there  is  only  one  way  of  invoking  them: 

Input  -  SUBROUTINE  INCTRL 

Analysis  -  SUBROUTINE  FLSOLV 
Rezoning  -  SUBROUTINE  RZMAIN 

The  remaining  two  modules  consist  of  utility  subroutines/  any  of 
which  may  be  called  as  needed: 

Output  -  SUBROUTINES  PRTXYZ ,  PRTVEL ,  PRZONE,  MSHMPO 
Restart  -  SUBROUTINES  FLSAVE ,  FLREST 
In  the  stand-alone  mode  of  operation/  a  driver  program  is 
needed  to  control  the  execution  of  the  input/  analysis/  output 
and  restart  modules.  The  rezoning  operation  is  normally  called 
as  needed  under  control  of  the  analysis  segment.  When  the  SBI 
analysis  is  coupled  with  a  structural  analysis  program/  the 
interface  subroutines  perform  this  control  function.  In  the 
coupled  mode/  the  structural  analyzer  acts  as  the  control 
program/  calling  the  SBI  routines  for  input/  output/  analysis  or 
restart  functions  as  needed. 

Communication  between  the  SBI  modules  and  a  calling  program 
is  accomplished  in  two  ways.  First/  a  real  array  is  supplied  to 
the  SBI  routines  as  working  storage,  which  is  dynamically 
allocated  as  needed.  Second,  a  series  of  ten  COMMON  blocks  is 
used  to  retain  key  solution  parameters.  The  contents  of  these 
blocks  are  described  below,  including  the  modules  in  which  they 
are  defined  and  used.  Single-letter  codes  used  to  identify  the 
modules  are:  (I)nput,  (A)nalysis,  (O)utput,  (R)estart,  re(Z)one. 
When  variables  in  COMMON  must  be  defined  in  the  driver  program  or 
the  structural  analyzer  interface,  a  'D'  code  appears. 
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COMMON  Block  /PROBID/ 


Name 


ITITLE 


COMMON 

Name 


NWORK 

LWORK 

LWA 

IFILL 

NUMNOD 

NUMEL 

NUMMAT 

NVAR 

NUMRWC 

NUMCON 

NREF 


NUMSBC 

NFILL 

LCORD 

LSVAR 

LCONN 

LMATL 

LRWBC 


Defined  Used  Description 


I  I AOR  Alphanumeric  title  ( CHARACTER*80 ) 


Block  /CTLBLK/ 

Defined  Used  Description 


D 

A 

I 

IZ 

IZ 

I 

I 


I 

I 

I 


I 

IZ 

IZ 

IZ 

I 

I 


I ARZ  Length  of  working  array  (real  words) 

IARZ  Length  of  working  array  (real  words) 

available  for  use  in  solution  phase 
IA  Last  word  allocated  during  input 

-  Reserved  space/  20  integer  words 

IAORZ  Current  number  of  nodes  in  mesh 
IAORZ  Current  number  of  elements  in  mesh 
IARZ  Number  of  materials  defined 

IAORZ  Number  of  solution  variables  (normally 

11-3  velocities/  density/  6  stresses/ 
and  effective  strain  parameter) 

IAR  Number  of  rigid-wall  boundary  conditions 
IAR  Number  of  contact  (moving  surface) 
boundary  conditions 

IZR  List  of  length  (3)  giving  the  number  of 
rezoning  reference  stations  in  each 
coordinate  direction. 

IAR  Number  of  symmetry  boundary  conditions 

—  Reserved  space/  29  integer  words 

IARZ  Location  of  coordinate  data  block  in  the 
main  working  array 

IARZ  Location  of  solution  variables  data  block 

IARZ  Location  of  element  connection  data  block 

IAR  Location  of  properties  data  block 

IAR  Location  of  rigid-wall  b.c.  data  block 
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Description 


Name  Defined  Used 


LCNBC 

I 

IAR 

Location 

Of 

contact  b.c.  data  block 

LBODY 

I 

IAR 

Location 

of 

body  forces  data  block 

LREFP 

I 

IAR 

Location 

of 

rezoning  reference  position 

data  block 

LSYBC 

I 

IAR 

Location 

of 

symmetry  b.c.  data  block 

LFILL 

- 

- 

Reserved 

space,  29  integer  words 

COMMON 

Name 

Block  /SOLPAR/ 

Defined  Used 

Description 

TIME 

IA 

AOR 

Current  time 

va  lue 

TMAX 

I 

AR 

Maximum  time 

value 

TREST 

IR 

R 

Last  time  at 

which  restart  file  created 

TPOST 

10 

OR 

Last  time  at 

which  results  file  created 

DTMIN 

I 

AR 

Lower  limit 

on  solution  time  step 

DTMAX 

I 

AR 

Upper  limit 

on  solution  time  step 

FILL 

- 

- 

Reserved  space/  10  real  words 

I  NCR 

IA 

AR 

Current  time 

step  number 

INCMAX 

I 

AR 

Upper  limit  < 

on  number  of  time  steps 

IREST 

I 

DR 

Flag  for  restart  output  (0=off/l=on) 

I  POST 

I 

DR 

Flag  for  results  output  (0=off/l=on) 

MFILL 

- 

- 

Reserved  space/  10  integer  words 
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COMMON  Block  /TOLBLK/ 


Name 

Defined 

Used 

Description 

DTFRAC 

I 

AR 

Fraction  of  critical  time  step  to  be 

VOLTOL 

I 

AR 

used  in  the  solution/  typically  0.7-0. 9 
Relative  volume  change  permitted  before 

4 

SDTTOL 

I 

AR 

rezoning  is  invoked,  typically  0.1-0. 5 
Relative  change  permitted  in  critical 

HGDAMP 

I 

AR 

time  step  before  rezoning  is  invoked, 
typically  0.25-0.5 

Hourglass  damping  fraction.  This  is  the 

DFACT 

I 

AR 

fraction  of  critical  damping  used  to 
suppress  hourglassing  instabilities, 
typically  0.05-0.5 

Time  step  reduction  due  to  hourglass 

damping.  DFACT  is  applied  to  the  stable 
time  step  to  preserve  stability.  Set  to 
SQRT ( 1 -HGDAMP* *2 ) -HGDAMP 


COMMON  Block  /FLSPAR/ 


Name 

Defined 

Used 

Description 

DT 

A 

AR 

Current  time 

step 

TQUIT 

D 

AR 

Maximum  time 

value  used  to  force  an  exit 

from  the  analysis  module  at  intermediate 
times  in  the  solution.  Overrides  TMAX . 

I  FIRST  D  AR  Flag  for  first  pass  into  analysis  module; 

forces  initial  conditions  to  be  applied. 
=1  for  first  pass/  =0  afterward 
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COMMON  Block  /FILNAM/ 


Name 

Defined 

Used 

Description 

NIN 

D 

IR 

Logical  unit 

number 

for 

data  input 

NOUT 

D 

IAORZ 

Logical  unit 

number 

for 

printing 

COMMON 

Block  /FLFILE/ 

Name 

Defined 

Used 

Description 

IUPOST 

D 

OR 

Logical  unit 

number 

for 

results  file 

IUREST 

D 

R 

Logical  unit 

number 

for 

restart  file 

IUWORK 

D 

RZ 

Logical  unit 

number 

for 

scratch  file  to 

be  used  during  rezoning.  The  file  must 
be  opened  for  unformatted  sequential 
access  prior  to  entering  either  the 
analysis  module  or  the  rezoning  module. 
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COMMON  Block  /FSINTR/ 


Name 

Defined 

Used 

Description 

NSNODE 

I 

AR 

Number  of  nodes  in  interface  set 

NSELEM 

I 

AR 

Number  of  interface  (contact)  elements 

NFSFIL 

D 

I  AR 

Direct  access  file  for  interface  data 

NIDREC 

D 

I  AR 

File  record 

for  interface  node  numbers 

NCNREC 

D 

I  AR 

File  record 

for  interface  element  data 

NNCREC 

D 

I  AR 

File  record 

for  interface  geometry  data 

NNVREC 

D 

I  AR 

File  record 

for  interface  velocity  data 

LSNID 

i 

AR 

Location  of 

interface  array  ISNID 

LSCON 

i 

AR 

Location  of 

interface  array  ISCON 

LXYZS 

i 

AR 

Location  of 

interface  array  XYZS 

LVELS 

i 

AR 

Location  of 

interface  array  VELS 

LTR 

i 

AR 

Location  of 

interface  array  TR 

LXMIN 

i 

AR 

Location  of 

interface  array  XMIN 

LXMAX 

i 

AR 

Location  of 

interface  array  XMAX 

LYMIN 

i 

AR 

Location  of 

interface  array  YMIN 

LYMAX 

i 

AR 

Location  of 

interface  array  YMAX 

LZMAX 

i 

AR 

Location  of 

interface  array  ZMAX 

LFSEND 

i 

AR 

End  of  data 

area  for  interface  arrays 

COMMON 

Name 

Block  /WRKSPC/ 

Defined  Used 

Description 

WORK  A  A  Scratch  space  for  fixed-dimension  arrays 

used  in  element  calculations.  Block 
must  be  at  least  136  words  in  length. 
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COMMON  Block  /PCOMl  / 


Name 

Defined 

Used 

Description 

PCOM 

A 

AZ 

Scratch  space  used  to  compress  rezoning 

station  data  prior  to  beginning  a 
rezone.  Block  must  be  at  least  180 
words  in  length. 


COMMON  Block  /DP  / 

Name  Defined  Used  Description 

IDP  D  I ARZ  Flag  for  DOUBLE  PRECISION  usage.  =1  if 

variables  of  type  REAL  are  of  the  same 
length  as  type  INTEGER;  =2  if  type  REAL 
variables  are  twice  as  long  as  INTEGERS. 


Typically,  parameters  NWORK ,  IFIRST,  NIN,  NOUT,  IUPOST, 
IUREST,  IUWORK  and  IDP  are  defined  in  the  driver  program  prior  to 
calling  any  of  the  SBI  segments.  The  intermediate  time  parameter 
TQUIT  is  defined  prior  to  entering  the  analysis  module  each  time 
it  is  called.  For  execution  in  the  coupled  mode,  the  following 
parameters  in  block  /FSINTR/  must  be  defined  prior  to  entering 
the  input  module:  NFSFIL,  NIDREC,  NCNREC ,  NNCREC ,  and  NNVREC . 

Entry  to  the  SBI  modules  requires  the  following  subroutine 
calls  (the  main  working  array  is  called  "A"): 
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1.  Input  Module 


CALL  INCTRL  (  A/  MODE  ) 

* 

where  MODE  =  0  for  stand-alone  execution#  1  for  coupled  mode 

2.  Analysis  Module 


CALL  FLSOLV  (  A ,  IFLAG,  IRZOPT,  MODE,  CD AT A  ) 

where  IFLAG  is  an  output  status  flag  (0  for  normal  completion, 
1  for  error(s)  in  element  geometry,  2  for  errors 
during  rezoning,  3  for  insufficient  storage) 
IRZOPT  >  0  forces  rezoning  upon  exit,  =  0  specifies 
rezoning  only  as  necessary 

MODE  =  0  for  stand-alone  execution,  1  for  coupled  mode 
CDATA  is  an  array  containing  interface  data  for  the 
coupled  mode  (see  subsection  4.8) 

3.  Output  Module 


CALL  PRTVEL  (  NUMNOD,  NVAR ,  A(LSVAR)  ) 

CALL  MSHMPO  (  NUMNOD,  A(LCORD),  NUMEL,  A(LCONN),  NVAR, 

A(LSVAR),  NREF ,  A(LREFP),  IUPOST,  INCR,  TIME  ) 

4.  Restart  Module 


CALL  FLSAVE  (  A,  IUREST  ) 
CALL  FLREST  (  A,  IUREST  ) 

5.  Rezoning  Module 


CALL  RZMAIN  (  NUMNOD,  NUMEL  ,  NVAR  ,  NREF  ,  A(LREFP), 

NWORK  ,  A  ,  LCORD ,  LSVAR,  LCONN  , 

I ERROR ,  IUWORK ,  NOUT  ) 

where  I  ERROR  is  an  output  status  flag  (see  subsection  4.6) 
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4.3  INPUT  MODULE 


The  input  module  of  the  SBI  code  controls  the  reading  of 
problem  data  describing  the  hydrodynamic  finite  element  mesh  (see 
Section  5.3).  During  the  input  procedure,  the  input  module 
performs  memory  allocation  for  data  as  it  is  read,  and  prestores 
addresses  at  which  the  mesh  interface  data  will  be  stored  in  the 
coupled  execution  mode. 

Entry  to  the  input  module  is  achieved  by  the  single 
subroutine  call 


CALL  INCTRL  (  A,  MODE  ) 


in  which  A  is  the  main  working  array,  and  MODE  is  a  flag 
indicating  the  execution  mode  (MODE=0  for  stand-alone,  >0  for 
coupled  mode) .  The  call  to  INCTRL  causes  the  following  input 
data  blocks  to  be  read  from  file  NIN: 


1.  TITL  -  problem  title 

2.  PARA  -  solution  parameters 

3.  TOLE  -  tolerances  for  rezoning  and  hourglass  forces 

4.  PROP  -  physical  properties 

5.  REFE  -  rezoning  reference  positions* 

6.  SYMM  -  symmetry  boundary  conditions 

7.  RIGI  -  rigid  wall  boundary  conditions 

8.  CONT  -  contact  (structure  interface)  data,  if  MODE>0 

9.  BODY  -  prescribed  body  forces 

10.  NODE  -  nodal  coordinate  data 

11.  INIT  -  initial  conditions 

12.  ELEM  -  finite  element  definitions 

Parameters,  tolerances,  sizing  data  and  address  information 
are  stored  in  COMMON  blocks  (see  subsection  4.2)  by  the  input 
module.  Remaining  data  are  stored  sequentially  in  the  main 
working  array  "A",  as  summarized  below. 


*  REFE  block  input  is  particular  to  rezoning  algorithm  I  (see 
subsection  2.6.2) ,  but  is  required  when  either  algorithm  is  used. 
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Location  Data  Type  Description 


LMATL 

Real 

Physical  properties  data 

LREFP 

Real 

Rezoning  reference  stations 

LSYBC 

Real 

Symmetry  boundary  condition  data 

LRWBC 

Real 

Rigid-wall  boundary  condition  data 

LBODY 

Real 

Body  forces 

LCORD 

Real 

Nodal  coordinates 

LSVAR 

Real 

Solution  variables,  with  initial  conditions 

applied 

LCONN 

Integer 

Element  connections 

Data  are  read  and  stored  in  the  above  sequence,  regardless  of  the 
order  of  data  blocks  in  the  input  file.  Note  that  data  blocks 
which  remain  unchanged  when  rezoning  is  performed  are  stored 
first.  The  coordinate,  solution  variable  and  element  connection 
blocks  may  change  in  size  when  a  rezone  occurs. 

In  the  coupled  mode,  the  input  module  also  processes 
interface  definition  data  ( CONT  data  block):  this  data  is 
organized  in  its  final  form  and  written  directly  to  random  access 
storage  for  use  by  the  interface  subroutines.  File  and  record 
numbers  are  defined  in  the  calling  program  and  are  stored  in 
COMMON  block  /FSINTR/ . 


4.4  RESTART  MODULE 


The  SBI  restart  module  is  used  to  save  or  restore  problem 
data  at  an  intermediate  stage  in  the  problem  solution.  Two 
subroutines  exist  for  this  purpose: 

-  FLSAVE  (save  data  for  restart)/  and 

-  FLREST  (restart  from  stored  data). 

Entry  to  these  routines  requires  the  following  subroutine  calls: 

-  CALL  FLSAVE  (  A,  LUNIT  ) 

-  CALL  FLREST  (  A,  LUNIT  ) 

Here  A  is  the  main  working  array  and  LUNIT  is  the  number  of  the 
sequential  file  to  which  the  restart  information  is  written. 

The  restart  file  consists  of  two  records/  containing  the 
contents  of  the  major  common  blocks  for  the  SBI  code  as  well  as 
the  main  working  array  A.  The  two  restart  records  are  written  as 
follows: 

WRITE  (LUNIT)  NWORK  ,  LWORK  ,  LWA  /  IFILL  / 


NUMNOD, 

NUMEL 

/ 

NUMMAT, 

NVAR 

/ 

NUMRWC 

NUMCON, 

NREF 

/ 

NUMSBC , 

NFILL 

/ 

LCORD  , 

LSVAR 

9 

LCONN  , 

LMATL 

9 

LRWBC 

LCNBC  , 

LBODY 

9 

LREFP  , 

LSYBC 

9 

LFILL 

TIME  , 

TMAX 

9 

TREST  , 

TPOST 

9 

DTMIN 

DTMAX  , 

FILL 

9 

INCR  , 

INCMAX 

9 

IREST  , 

I  POST 

9 

MFILL 

DTFRAC , 

VOLTOL 

9 

SDTTOL , 

HGDAMP 

9 

DFACT 

DT  , 

TQUIT 

9 

IFIRST 

LASTWD  =  LCONN  +  8*NUMEL/IDP  -  1 
WRITE  (LUNIT)  (A(l),I=l , LASTWD ) 

Interface  data  for  operation  in  the  coupled  mode  must  be  saved 
elsewhere  for  restart,  since  the  data  required  for  linking  the 
SBI  solution  with  a  particular  structural  analyzer  will  be  highly 
program-dependent. 
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4.5  SOLUTION  MODULE 


The  SBI  solution  module  performs  a  segment  of  the  explicit/ 
hydrodynamic  solution,  consisting  of  a  specified  time  interval  or 
a  predetermined  number  of  time  steps  (whichever  comes  first). 

The  organization  of  input  and  output  data  is  as  described  in 
subsections  4.3  and  4.7.  However,  since  the  solution  package  can 
invoke  the  rezoning  procedure,  the  input  and  output  array  sizes 
may  differ. 

The  solution  module  is  entered  through  the  subroutine  call: 

CALL  FLSOLV  (  A,  I FLAG ,  IRZOPT,  MODE,  CDATA  ) 

Formal  parameters  for  the  module  are  as  follows: 

A  =  Main  working  array,  described  in  subsection  4.3. 

IFLAG  =  Output  status  flag 

=0,  Normal  completion  ( TQUIT  or  INCMAX  reached) 

=1,  Error(s)  in  element  geometry  detected 

=2,  Error(s)  occurred  during  rezoning  operation 

=3,  Working  storage  allocation  exceeded 

IRZOPT  =  Switch  for  rezoning  of  final  solution 

=0,  Rezoning  is  performed  only  as  necessary 

>0,  Forces  rezoning  before  leaving  solution  module 

MODE  =  Switch  for  execution  mode 
=0,  Stand-alone 

=1,  Coupled  with  structural  analysis 

CDATA  =  Structural  interface  data,  used  only  when  MODE=l 

The  actual  problem  data  are  supplied  to  the  solution  routines  in 
the  work  array  A,  and  in  the  COMMON  blocks  described  in 
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subsection  4.2.  The  contents  of  array  CDATA,  which  describe  the 
structure  model  interface/  are  summarized  in  subsection  4.8. 

A  hydrodynamic  solution  may  be  performed  either  through  a 
single  call  to  FLSOLV /  or  by  a  series  of  calls;  minor  differences 
may  occur  between  these  two  methods  when  rezoning  is  performed 
(particularly  if  IRZOPT>0) .  The  usual  reasons  for  performing  a 
series  of  calls  to  the  solution  module  are  to  allow  results  and 
restart  data  to  be  stored  at  intervals  throughout  the  analysis/ 
or  to  force  results  to  be  generated  at  regular  intervals  during 
the  solution. 
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4.6  REZONING  MODULE 

The  SBI  rezoning  module  performs  a  complete  reorganization 
of  the  finite  element  grid  for  the  hydrodynamic  solution.  This 
module  is  suitable  for  use  with  other  forms  of  finite  element 
data  as  well,  provided  it  conforms  to  the  input  conventions 
described  below.  In  effect,  the  input  to  the  rezoner  is  a  finite 
element  model  composed  of  eight-node  hexahedral  elements, 
together  with  certain  information  to  control  the  rezoning 
process;  the  output  is  a  different  finite  element  model  with 
minimally-distorted  elements. 

There  are  two  primary  limitations  which  exist  in  the 
rezone  module.  First,  the  model  must  be  composed  entirely  of 
eight-node  brick  elements.  Second,  the  number  of  nodal  solution 
variables  is  limited  to  11  by  fixed  DIMENSION  statements  in 
several  of  the  rezoning  subroutines. 


Entry  to  the  rezoning  package  using  rezoning  algorithm  I 
(subsection  2.6.2)  requires  a  single  subroutine  call  having  the 
form: 


CALL  RZMAIN  (  NUMNOD,  NUMEL  ,  NVAR  ,  NREF  ,  XREF  , 
NWORK  ,  A  ,  LCORD  ,  LSVAR  ,  LCONN  , 

IERROR,  LUNIT  ,  NOUT  ) 

Formal  parameters  for  the  routine  are: 


NUMNOD  (  input) = 
(output)  = 
NUMEL  (  input) = 
(output)  = 
NVAR  (  input) = 


Number  of  nodes  in  input  mesh 
Number  of  nodes  in  output  mesh 
Number  of  elements  in  input  mesh 
Number  of  elements  in  output  mesh 
Number  of  solution  variables  per  node. 


0  <  NVAR  <12 


NREF  (  input) =  Integer  array  of  length  3  containing  the 

number  of  reference  positions  in  the  Xl, 

X2,  and  X3  directions,  respectively 
XREF  (  input) =  Real  array  containing  the  reference  station 

positions  in  ascending  order  along  each  of 
the  coordinate  directions.  XREF(i,j)  is 
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NWORK 

A 


LCORD 

LSVAR 

LCONN 

IERROR 


the  j-th  reference  station  in  the  i-th 
direction.  The  leading  dimension  of  XREF 
must  be  MAX(NREF(1)  ,NREF(2)  ,NREF(3)  ) 

(  input) =  Number  of  floating-point  words  available  in 
the  working  array  'A'  (see  below) 

(  input) =  Primary  data  array,  containing  the  nodal 
coordinates,  solution  variables,  and 
element  connections,  at  the  starting 
locations  described  below  (see  LCORD, 

LSVAR,  and  LCONN) .  These  arrays  must  be 
stored  in  the  order  stated  above. 

(output) =  Primary  data  array,  containing  the  same 
data  for  the  new  finite  element  mesh. 

(  input)3  Starting  location  for  nodal  coordinate  data 
in  array  'A'  for  the  input  mesh 
(output)3  Starting  location  for  nodal  coordinate  data 
in  array  'A'  for  the  output  mesh 
(  input) 3  Starting  location  for  solution  variables 
data  in  array  'A'  for  the  input  mesh 
(output) 3  Starting  location  for  solution  variables 
data  in  array  'A'  for  the  output  mesh 
(  input) 3  Starting  location  for  element  connection 
data  in  array  'A'  for  the  input  mesh 
(output)3  Starting  location  for  element  connection 
data  in  array  'A'  for  the  output  mesh 
(output) 3  Exit  error  code,  defined  as  follows: 

3  0  ,  Normal  completion 

3  1  ,  Number  of  nodes  illegal 

3  2  ,  Number  of  elements  illegal 

3  3  ,  Number  of  solution  variables  illegal 

3  4  ,  Number  of  reference  stations  illegal 

for  one  or  more  coordinate  directions 
3  5  ,  Reference  station  values  for  one  or 

more  directions  not  in  ascending 
order 

3  10  ,  Input  data  (LCORD,  LSVAR,  LCONN)  out 

of  sequence 
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=  20  ,  Parameters  LCORD,  LSVAR  and/or  LCONN 

inconsistent  with  number  of  nodes, 
elements  or  solution  variables 
declared 

=  100  ,  Storage  exceeded  before  rezoning 

=  1000,  Storage  exceeded  during  generation  of 
new  node  points 

=  1010,  Storage  exceeded  during  generation  of 
new  elements 

=  2000,  Error (s)  detected  during  element  fill 
operation 

LUNIT  (  input) =  Logical  unit  number  of  scratch  file  for  use 

in  the  rezoner.  This  file  must  be  opened 
for  unformatted,  sequential  I/O  prior  to 
calling  RZMAIN,  and  may  be  discarded  or 
overwritten  afterward. 

NOUT  (  input) =  Logical  unit  number  for  line  printer 

output,  including  error  messages  and 
rezoning  statistics.  If  NOUT=0,  no  printed 
output  is  generated  by  the  rezoning  module. 

No  information  is  passed  to  the  rezoning  package  through  COMMON, 
to  facilitate  its  use  as  a  utility  module  in  other  applications. 

Using  rezoning  algorithm  II,  entry  to  the  rezoning  module 
requires  the  following  subroutine  call: 


CALL  RZMAIN  (NUMNOD,  NUMEL  ,  NPASS,  LENGTH,  IERROR, 
NVAIR  ,  XYZ  ,  VAR  ,  NCON  ,  WORK  , 

IWORKl ,  I WORK 2 ,  LUNIT) 


Formal  parameters  for  the  routine  are: 


NUMNOD 

(input) 

=  Number  of 

nodes  in 

mesh 

NUMEL 

( input) 

=  Number  of 

elements 

in  mesh 

NPASS 

(input) 

=  Number  of 

rezoning 

passes  to 

be  performed 

(normally 

1) 

LENGTH 

( input) 

=  Length  of 

work  array  IWORK2 

(see  below 
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IERROR 


(output) 


NVAR 

( input) 

XYZ 

( input) 

(output) 

VAR 

( input) 

(output) 

NCON 

( input) 

WORK 

( input) 

IWORK1 

( input) 

IWORK2 

( input) 

LUNIT 

( input) 

Exit  error  code,  defined  as  follows: 

=0 ,  normal  completion 
=1,  error (s)  in  remapping  solution 
variables 

-I#  length  of  IWORK2  must  be  increased  by 
n  words  to  perform  rezoning 
Number  of  solution  variables  (normally  11) 
Cartesian  coordinates  of  nodes  (3xNUMNOD) 
Coordinates  in  rezoned  mesh 
Solution  variables  array  (NVAR  x  NUMNOD) 
Remapped  solution  variables 
Element  connectivity  (8  x  NUMEL) 

Work  array  at  least  3  x  NUMNOD  words  in 
length 

Work  array  at  least  NUMNOD  integer  words 
in  length 

Work  array  of  estimated  length  at  least  8 
x  NUMEL  integer  words 
Logical  unit  number  of  scratch  file 
available  for  use  by  RZMAIN.  This  file 
must  be  opened  for  unformatted,  sequential 
I/O  prior  to  calling  RZMAIN. 
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4.7  OUTPUT  MODULE 


The  SBI  output  module  consists  of  four  subroutines  which  may 
be  called  at  any  time  to  produce  printed  or  file  output.  Three 
of  these  (PRTXYZ,  PRTVEL ,  MSHMPO)  are  used  routinely  for  output 
of  analysis  results;  the  fourth,  PRZONE,  is  a  diagnostic  routine 
which  summarizes  intermediate  results  during  rezoning. 

PRTXYZ  prints  a  listing  of  current  nodal  coordinates  on  a 
designated  output  file.  Since  the  hydrodynamic  mesh  is  updated 
at  each  time  step,  this  output  describes  the  deformed  geometry  of 
the  finite  element  model  at  the  current  time.  PRTXYZ  is  called 
as  follows: 

CALL  PRTXYZ  (  NUMNOD,  XYZ ,  NOUT  ) 

in  which  NUMNOD  is  the  number  of  nodes  currently  defined  for  the 
mesh,  XYZ  is  the  coordinate  data  in  a  3  by  NUMNOD  array,  and  NOUT 
is  the  logical  unit  number  for  output. 

PRTVEL  is  similar  to  PRTXYZ,  but  prints  velocity  information 
in  a  similar  format.  The  only  important  differences  are  in  the 
header  information  printed  by  the  routines,  and  in  the  storage 
format  expected.  Because  the  velocities  are  normally  stored 
together  with  other  solution  variables,  the  calling  sequence  for 
PRTVEL  includes  NVAR,  the  number  of  solution  variables;  the  first 
three  entries  of  each  row  of  the  NVAR  by  NUMNOD  array  VAR.  The 
call  to  PRTVEL  looks  like: 

CALL  PRTVEL  (  NUMNOD,  NVAR,  VAR  ) 

Output  is  written  to  the  unit  number  NOUT  as  defined  in  COMMON 
block  /FILNAM/. 

MSHMPO  generates  postprocessing  file  output  in  a  form 
similar  to  the  MAGNA  results  file  "MPOST"  [16].  The  calling 
sequence  for  MSHMPO  is: 

CALL  MSHMPO  (  NUMNOD,  XYZ,  NUMEL,  NCON  ,  NVAR,  VAR, 

NREFX  ,  XR  ,  NOUT  ,  INCRNO,  TIME  ) 

Parameters  NUMNOD,  XYZ,  NVAR,  VAR  and  NOUT  are  as  described  above 
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for  the  other  output  routines.  NUMEL  defines  the  number  of 
elements/  and  NCON (8 /NUMEL )  contains  the  element  connections. 
NREFX(3)  contains  the  number  of  rezoning  reference  stations  in 
each  coordinate  direction,  and  XR  the  actual  reference  station 
coordinates.  The  minimum  and  maximum  reference  position  values 
are  written  as  eight  extra  nodes  with  the  postprocessing  file 
output,  so  that  plotting  routines  which  examine  the  coordinates 
for  minimum  and  maximum  coordinate  values  to  determine  scaling  of 
plots  will  produce  sequences  of  plots  having  the  same  scaling 
characteristics.  INCRNO  and  TIME  define  a  time  step  number  and 
time  value  to  be  included  in  the  output  headers  on  the  results 
file. 

Subroutine  PRZONE  is  normally  used  only  for  debugging.  It 
must  be  called  from  subroutine  REZONE  in  the  rezoning  module,  by 
modifying  the  source  code  in  that  module. 
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4.8  MESH  INTERFACE  CODE 


The  mesh  interface  utilities  consist  of  several  subroutines 
which  link  the  SBI  code  to  a  nonlinear  structural  analysis 
program,  for  the  purpose  of  performing  hydrodynamic  impact 
analysis.  These  routines  fall  into  three  categories: 

(a)  Interface  data  input; 

(b)  Interface  geometry  construction;  and 

(c)  Analysis  interface. 


Each  group  of  utility  routines 
descriptions  are  intended  only 
the  SBI  code  with  a  structural 
details  of  doing  so  are  almost 


are  described  briefly  below.  The 
as  a  general  guide  to  interfacing 
analysis  program,  since  the 
completely  program-dependent. 


The  interface  data  input  routines  include  FSDATA ,  FSINPT, 

FSICON ,  and  FSNIDS.  FSDATA  acts  as  a  driver  for  the  SBI  input 
module,  opening  the  necessary  files  and  defining  necessary 

parameters  in  COMMON  (see  subsection  4.2).  The  remaining 
routines  are  called  from  the  input  module  (subsection  4.3) 
whenever  MODE=l,  to  read  the  structural  interface  data.  The  data 
defined  in  these  input  routines  consists  of  interface  element 
definitions  in  terms  of  the  nodes  of  the  structural  model  (see 
subsections  2.7  and  5.3  for  details).  The  structural  analysis 
program  typically  calls  FSDATA  once  during  the  input  phase  of 
execution. 


Interface  geometry  construction  is  performed  in  subroutines 
FSPUTC ,  FSGETC ,  and  FSGETU .  FSPUTC  is  called  by  the  structures 
program  to  save  the  nodal  coordinates  on  file  as  soon  as  they  are 
available.  FSGETC  must  be  called  later,  after  the  interface  data 
have  been  read.  The  result  of  this  stage,  which  is  performed 
only  once,  is  to  save  original  coordinates  for  the  interface 
elements.  These  will  be  used  repeatedly  in  the  analysis  phase  to 
compute  current  positions  of  the  interface  elements.  This  phase 
is  unnecessary  if  a  relative  (updated  Lagrangian)  description  of 
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motion  is  used  in  the  structural  analysis,  since  updated  nodal 
coordinates  should  be  available  at  any  time  during  the  solution. 

The  analysis  interface  consists  of  several  subroutines 
(FLCALC,  FSXYZV ,  FSTRAN ,  FSCONT,  and  FSFORC)  which  are  called  at 
each  major  (structure)  time  step  of  the  coupled  solution.  The 
primary  entry  point  is  FLCALC,  which  controls  the  assembly  of 
structure  mesh  data  for  the  impact  analysis,  calls  the  solution 
module,  and  controls  the  retrieval  of  interface  nodal  forces. 

The  secondary  routines  are  described  briefly  below. 


FSXYZV 

FSTRAN 


FSCONT 


FSFORC 


retrieves  coordinate  and  velocity  information  for 
structure  nodes  which  lie  on  the  interface, 
performs  preliminary  geometric  calculations  for  all 
interface  elements  based  on  their  current  geometry; 
these  include  transformations  to  local  coordinates, 
and  storing  of  minimum  and  maximum  coordinate 
values. 

this  routine  is  called  from  FLSOLV  (solution 
module)  whenever  MODE>0,  to  perform  the  actual 
interface  contact  calculations.  The  results  are 
coordinates  and  velocities  updated  to  reflect  the 
existing  contact  conditions,  and  nodal  values  of 
the  interface  forces. 

accepts  the  interface  forces  computed  in  FSCONT, 
and  generates  nodal  forces  to  be  applied  to  the 
structural  mesh  during  the  current  time  step. 


The  analysis  interface  routines  FSDATA  and  FLCALC  also  use  the 
restart  module  subroutines  to  save  and  restore  hydrodynamic  mesh 
data  as  needed  during  the  solution. 


Interface  data  which  is  passed  to  the  solution  module  during 
a  coupled  analysis  is  arranged  sequentially  in  the  array  CDATA 
(see  subsection  4.5).  This  data  is  assembled  by  FLCALC  (above), 
and  contains  the  following  lists  and  arrays: 
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ISNID 


Structure  model  nodes  by  sequential  order  as  they 
appear  on  the  interface 

ISCON  -  Interface  panel  connections  (4  per  element) 

XYZS  -  Interface  nodal  coordinates 

VELS  -  Interface  nodal  velocities 

TR  -  Global-to-local  transformation  for  each  interface 

element 

XMIN  -  Minimum  local  'X'  coordinate  for  each  interface 
element 

XMAX  -  Maximum  local  'X'  coordinate  for  each  interface 
element 

YMIN /  YMAX ,  ZMAX  -  (Similar  to  XMIN  and  XMAX) 

H FORCE  -  Area  for  output  of  interface  nodal  forces 

All  of  the  above  arrays  except  H FORCE  are  input  to  the  solution 
module:  there  size  is  solely  dependent  upon  the  number  of  nodes 
and  elements  on  the  structural/hydrodynamic  interface,  and 
therefore  does  not  change  over  the  course  of  the  solution.  The 
interface  element  data  is  prestored  before  the  contact  analysis 
to  save  time  recomputing  these  parameters,  at  the  expense  of  some 
array  storage  space.  The  total  array  storage  necessary  is 
10x(number  of  nodes  on  interface)  +  18x(number  of  interface 
elements) . 
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SECTION  5 
PROGRAM  OPERATION 

This  section  provides  a  brief  description  of  procedures  for 
operating  the  soft  body  impact  analysis  code.  Both  stand-alone 
and  coupled  modes  of  execution  are  discussed;  the  coupled  mode 
is  described  with  reference  to  a  particular  structural  dynamics 
program#  MAGNA/  which  has  been  linked  with  the  hydrodynamic 
analysis  during  the  present  project.  Subsection  5.3  summarizes 
the  input  blocks  and  parameters  necessary  for  execution  in  either 
mode. 
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5.1  OPERATION  IN  STAND-ALONE  MODE 

For  stand-alone  execution  of  the  hydrodynamic  analysis,  it 
is  necessary  only  to  provide  the  proper  input  files  to  the  code 
and  to  initiate  execution.  Restart  and  postprocessing  files  must 
be  saved  on  some  computers  after  execution  is  completed.  Some 
examples  of  job  control  language  for  CDC,  VAX  and  CRAY  computer 
systems  are  listed  below. 

Example  1.  CDC-6600  or  CYBER-175  under  NOS  (SUBMIT  file) 

/JOB 

jobname. 

/USER 

CHARGE,*. 

GET , INDAT= filename. 

GET , SBILG0/UN=D8  20139 . 

DEFINE, MPOST. 

SETTL , nnnnn. 

SBILGO. 

/EOR 

/EOF 

Example  2.  VAX-11/780  under  VAX/VMS  ( DCL  command  file) 

$  SET  VERIFY 

$  ASSIGN /US ER_MODE  input_f ile_spec  FOR005 

$  ASSIGN/USER _ MODE  output_f i le_spec  FOR006 

$  ASSIGN/USER_MODE  post_file  spec  MPOST 
$  RUN  SBI 
$  LOGOUT 
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Example  3.  CRAY-l/S  or  CRAY-X/MP  under  COS 


JOB,JN= j obname , Tnnnn . 

ACCOUNT , AC=a  cc t no , US= u se rn o , UPW  =  pas sword. 

ACCESS, DN=INDAT,PDN=f ilenamel. 

ASSIGN , A=FT9  5 , DN=MPOST . 

ACCESS, DN=S BI ,PDN=SBI . 

SBI. 

SAVE,DN=MPOST,PDN=f ilename2. 

EOF 

Example  4.  CRAY-l/S  under  COS  (United  Information  Services) 

JOB,Tnnnn. 

ACCOUNT , userno, password. 

ASSIGN , A=FT9  5 , DN=MPOST . 

GET, INDAT=f ilenamel /CI=TTY. 

GET, SBI . 

SBI . 

PUT , MPOST/CO=BAT/D . 

DFD , SBI DAY , R . 

PUT , $OUT=SBIOUT/CO=BAT . 

EOF 

In  the  above  job  control  files,  the  file  INDAT  contains 
problem  input  data  in  the  form  described  in  subsection  5.3.  The 
results  file  is  called  MPOST  in  all  of  the  examples.  Since  our 
primary  interest  is  in  the  coupled  mode  execution,  no  standard 
procedures  for  accessing  restart  files  have  been  developed.  The 
driver  for  the  hydrodynamics  program  can  be  modified  to  include 
the  restart  function;  for  each  example  above,  appropriate  ASSIGN, 
GET,  SAVE  or  other  file  control  statements  must  be  added  to 
define  and  save  the  restart  file. 
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5.2  OPERATION  IN  COUPLED  MODE 

When  the  soft  body  impact  analysis  is  coupled  with  a  finite 
element  structural  model,  three  additional  steps  are  necessary: 

(1)  provide  the  input  file  for  the  hydrodynamic  mesh; 

(2)  activate  the  appropriate  option  in  the  structures 
program  input;  and 

(3)  add  the  necessary  file  control  statements  to  the  job 
control  file. 

In  the  case  of  coupled  execution  with  MAGNA,  the  input  file 
for  the  soft  body  mesh  (subsection  5.3)  must  be  provided  as  the 
formatted  file  "FSDATA" .  The  input  data  must  contain  the  CONT 
input  block,  which  defines  the  structural/hydrodynamic  mesh 
interface . 

The  option  flag  in  MAGNA  which  activates  the  coupled  impact 
analysis  is  I0PT(18).  IOPT  option  flags  are  discussed  in  Section 

8.2  of  the  MAGNA  user's  manual.  Necessary  values  of  this  and 
other  IOPT  parameters  are  listed  below. 

IOPT  Columns  Value  Description 


2  5-8  2  Selects  transient  dynamic  solution 

3,4  9  -16  2  Selects  nonlinear  analysis 

18  68  -72  1  Selects  coupled  impact  analysis 

All  other  MAGNA  options  and  parameters  may  be  switched  on  or  off 
as  needed  to  define  the  particular  problem  at  hand. 

The  only  additional  job  control  statement  needed  for  the 
coupled  analysis  is  a  single  statement  to  access  the  data  file 
FSDATA  prior  to  execution;  for  example: 

GET, FSDATA= filename .  (CDC  under  NOS) 

$  ASS IGN/USER_MODE  file_spec  FSDATA  (VAX  under  VMS) 
ACCESS , DN= FS  DATA, PDN=  f i lename . 
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ASSIGN, A=FT40,DN=FSDATA. 


(CRAY-1  under  COS) 


GET, FSDATA=f ilename/CI=TTY. 

ASSIGN, A=FT40,DN=FSDATA.  (CRAY-1  at  U.I.S.) 

The  restart  and  postprocessing  files  generated  in  the  soft  body 
analysis  are  merged  automatically  with  the  corresponding  MAGNA 
files  (NOREST,  NRSTAP ,  and  MPOST)  . 

In  coupled  execution  mode,  additional  data  are  also  output 
to  the  MPOST  (postprocessor)  file  to  provide  a  less  complicated 
summary  of  the  hydrodynamic  solution  and  to  save  the  computed 
values  of  interface  forces.  These  data  are  interspersed  with  the 
data  blocks  which  usually  appear  on  the  MPOST  file,  and  have  the 
following  identifying  keywords: 

-  CMOM  :  center-of -mass  and  momentum  data  for  the 

hydrodynamic  mesh; 

-  HYDF  :  hydrodynamic  impact  forces  by  structure  node 

number;  and 

-  TOTF  :  total  hydrodynamic  forces  in  all  directions. 

The  CMOM  data  is  useful  for  visualizing  the  gross  motion  of  the 
soft  body  mesh,  and  for  evaluating  the  momentum  transfer  to  the 
target.  The  HYDF  data  blocks  provide  sufficient  information  to 
reconstruct  the  force-versus- time  history  generated  in  a  fully 
coupled  solution,  so  that  later  simulations  with  minor  parameter 
changes  can  be  run  in  an  uncoupled  mode  using  specified  force 
data.  The  force  resultants  from  the  TOTF  blocks  are  useful  when 
compared  with  momentum  data  for  the  hydrodynamic  model,  since 
the  combination  of  the  two  provides  a  means  of  verifying  that  the 
overall  motion  of  the  soft  body  is  correct. 

Output  data  and  formats  for  the  additional  data  blocks 
mentioned  above  are  as  follows: 
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Block 

Line  (s ) 

Data  Description 

Format 

CMOM 

1 

Block  header 

A4f  IX 

SBI  Increment 

no . 

15 

Time 

El5 . 7 

X  position  of 

C.M. 

E15.7 

Y  position  of 

C.M. 

El 5. 7 

Z  position  of 

C.M. 

El 5. 7 

X  direction  momentum 

El  5 . 7 

Y  direction  momentum 

E15.7 

Z  direction  momentum 

El 5. 7 

Block 

Line (s ) 

Data  Description 

Format 

HYDF 

1 

Block  header 

A4,  IX 

Structure  incr .  no. 

15 

Time 

E15.7 

Structure  node  no. 

15 

X  direction  force 

El  5.7 

Y  direction  force 

El 5. 7 

Z  direction  force 

E15.7 

HYDF 

2,3,... 

Blank 

2  5X 

Structure  node  no. 

15 

X  direction  force 

El 5. 7 

Y  direction  force 

El 5. 7 

Z  direction  force 

El 5. 7 
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Block 

Line (s ) 

Data  Description 

Format 

TOTF 

1 

Block  header 

A4 ,  IX 

Structure  incr.  no. 

15 

Time 

El  5 . 7 

Blank 

5X 

X  force  (total) 

El 5. 7 

Y  force  (total) 

El  5.  7 

Z  force  (total) 

El 5. 7 

The  MPOST  data  blocks  described  in  this  Section  can  be  extracted 
from  the  complete  postprocessing  file  simply  by  searching  for  the 
appropriate  block  headers.  It  should  be  noted  that  the  HYDF  data 
block  varies  in  length,  depending  on  the  number  of  structure 
nodes  which  experience  non- zero  impact  forces.  The  end  of  the 
HYDF  block  is  signalled  only  by  the  beginning  of  the  TOTF  block, 
which  consists  of  a  single  line.  If  no  impact  forces  exist  for  a 
particular  time  increment,  only  the  TOTF  line  will  appear  in  the 
file. 
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5.3  INPUT  DATA  DESCRIPTIONS 


Problem  input  for  the  hydrodynamic  analysis  is  supplied  on  a 
single  formatted,  sequential  file.  Data  are  separated  into  input 
"blocks";  each  input  block  begins  with  a  four  (or  more)  character 
identifier,  followed  by  the  appropriate  data;  most  end  with  a 
single  blank  line.  Input  blocks  may  appear  on  the  file  in  any 
order.  Comment  lines  may  be  inserted  in  the  input  file  after  the 
terminator  (blank  line)  for  any  block,  or  at  the  beginning  of  the 
file.  Comment  lines  have  no  special  format,  but  should  begin 
with  a  character  other  than  A-Z  to  avoid  duplicating  valid  input 
block  names. 


Data  blocks  recognized  by  the  input  module  are  as  follows; 


TITL  -  Problem  title 

NODE  -  Nodal  coordinate  data 

ELEM  -  Finite  element  connection  data 

PROP  -  Physical  properties  data 

BODY  -  Body  forces 

INIT  -  Initial  conditions 

SYMM  -  Symmetry  or  antisymmetry  boundary  conditions 
RIGI  -  Rigid  wall  boundary  conditions 
PARA  -  Solution  parameters  and  options 
TOLE  -  Tolerance  parameters 

REFE  -  Reference  position  data  for  rezoning 

CONT  —  Contact  interface  data  (coupled  execution  only). 


The  content  and  format  of  each  of  these  input  blocks  are 
described  in  the  remainder  of  this  subsection. 


In  general,  the  names  of  input  data  items  follow  FORTRAN 
typing  conventions;  that  is,  names  beginning  with  letters  I-N  are 
of  integer  data  type,  and  names  beginning  with  A-H  or  0-Z  are 
real  (floating-point)  data.  Where  character  (alphanumeric)  input 
is  required  (as  in  the  block  headers  and  title  information),  this 
is  noted  in  the  data  description.  Integer  data  may  be  entered 
anywhere  within  the  fields  provided,  but  may  not  contain  decimal 
points  or  exponents  (E).  Floating-point  data  may  be  placed 


77 


anywhere  within  the  data  field/  and  may  be  expressed  in  virtually 
any  format.  For  example,  20,  20.0,  2E1,  an  .2E2  all  represent  a 
numeric  value  of  "20". 
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BODY  FORCES  Input  Block 


Header  :  BODY 
End  :  (1  Record) 

Record  Format  #1  : 

\ 

.  .  .  .10  .  .  .  .20  .  .  .  .30 

XBODY  YBODY  ZBODY 

Input  Variables  : 


XBODY 

=  X 

direction 

body 

force 

YBODY 

=  Y 

direction 

body 

force 

ZBODY 

=  Z 

direction 

body 

force 

Notes  : 

(1)  -  Body  forces  are  defined  in  global  Cartesian  coordinate 

directions,  and  are  constant  throughout  the  solution. 

(2)  -  Body  force  data  is  expressed  in  terms  of  force  per  unit 

mass.  For  example,  gravity  acting  in  the  negative  "Y" 
direction  is  specified  by  setting  YBODY  =  -g,  where  g 
is  the  acceleration  of  gravity  in  units  consistent  with 
the  other  problem  data. 
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CONTACT 


INTERFACE  Input  Block 


Header  :  CONT 

End  :  blank  line 

Record  Format  #1  : 

.  .  .  .  5 

NSELEM 

Record  Format  #2  :  (repeat  as  needed  to  define  NSELEM  interface 

elements) 

.  .  5  .  .10  .  .15  .  .20  .  .25  .  .30  .  .35 

I EL  I EG  ING  N ( 1 )  N(2)  N(3)  N(4) 


Input  Variables  : 


NSELEM 
I  EL 
I  EG 
ING 
N(i) 


Total  number  of  interface  elements  to  be  defined 
Interface  element  number  ( 1 , 2 ,  . . . , NSELEM) 

Element  number  increment  for  use  in  generation 

Node  number  increment  for  use  in  generation 

i-th  connected  node  for  interface  element  'IEL't  note 

that  node  numbers  refer  to  the  nodes  of  the  structural 

model 


Notes  : 

(1)  -  This  input  block  defines  that  portion  of  the  surface  of 
the  structural  finite  element  model  with  which  the  soft 
body  (hydrodynamic)  model  is  expected  to  come  into 
contact. 
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Interface  elements  normally  correspond  to  exterior  faces 
of  solid  or  shell  elements  in  the  structures  model. 

Node  numbers  which  define  the  interface  elements  are 
specified  in  terms  of  the  node  numbers  in  the  structure 
finite  element  model/  not  the  hydrodynamic  model. 

Each  interface  element  is  defined  by  four  corner  nodes/ 
which  appear  in  counter-clockwise  order  when  viewed 
from  the  outside  of  the  contact  surface. 
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ELEMENTS 


Input  Block 


Header  :  ELEM 

End  :  Blank  line 

Record  Format  #1  : 

.  .  5 

NUMEL 

Record  Format  #2  :  (repeat  as  needed  to  define  NUMEL  elements) 

•  •  5  .  .10  •  .15  •  .20  .  .25  •  .30  •  .35  •  .40  •  .45  •  .50  •  .55 

IEL  I EG  ING  N ( 1 )  N ( 2 )  N(3)  N(4)  N(5)  N(6)  N(7)  N(8) 


Input  Variables  : 


NUMEL 
IEL 
I  EG 
ING 
N(i) 


Total  number  of  elements  to  be  defined 
Element  number  (1/2/..., NUMEL ) 

Element  number  increment  for  use  in  generation 
Node  number  increment  for  use  in  generation 
i-th  connected  node  for  element  'IEL' 


Notes  : 


(1)  -  The  relative  positions  of  local  nodes  1-8  of  an  element 
are  indicated  in  the  diagram  below. 


(4) - (3) 


83 


NODES 


Input  Block 


Header  :  NODE 

End  :  Blank  line 

Record  Format  #1  : 


...  5 

NUMNOD 

Record  Format  #2  :  (repeat  as  needed  to  define  NUMNOD  nodes) 

.  .  .  5  .  .  .10  . 20  . 30  . 40 

NODE  INGEN  XCORD  YCORD  ZCORD 

Input  Variables  : 

NUMNOD  =  Total  number  of  nodes  to  be  defined 
NODE  =  Node  number 

INGEN  =  Node  number  increment  for  use  in  generation 
XCORD  =  X  coordinate 
YCORD  =  Y  coordinate 
ZCORD  =  Z  coordinate 

Notes  : 

(1)  -  Valid  node  numbers  are  integers  from  1  through  NUMNOD. 

(2)  -  Not  all  nodes  between  1  and  NUMNOD  need  be  defined: 

however,  storage  is  allocated  for  NUMNOD  nodes,  and  the 
omission  of  large  numbers  of  nodes  may  lead  to  excessive 
storage  requirements. 

(3)  —  Nodal  data  may  be  read  or  generated  in  any  order. 
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initial 


CONDITIONS  Input  Block 


* 


f 


Header  :  INIT 
End  :  (l  Record) 

Record  Format  #1  : 

.  .  .  .10  .  .  .  .20  .  .  .  .30  .  .  .  .40 

XVEL  YVEL  ZVEL  DENS 


Input  Variables  : 


XVEL 

YVEL 

ZVEL 

DENS 


=  Initial 
=  Initial 
=  Initial 
=  Initial 


veloci ty 
ve loci ty 
velocity 
value  of 


in  X  direction 
in  Y  direction 
in  Z  direction 
material  density 


Notes  : 

(1)  -  Initial  conditions  specified  in  this  block  apply  to  all 

nodes  of  the  model, 

(2)  -  Except  in  special  circumstances,  the  material  density 

appearing  above  should  agree  with  the  reference  density 
defined  in  the  PROPerty  input  block. 
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PARAMETERS 


Input  Block 


Header  :  PARA 

End  :  (2  Records) 

Record  Format  #1  : 

.  .  .  .10  .  .  .  .20  .  .  .  .30  .  .  .  .40  .  .  .  .50  .  .  .  .60 

TIME  TMAX  TREST  TPOST  DTMIN  DTMAX 

Record  Format  #2  : 

.  .  .  5  .  .  .10  .  .  .15  .  .  .20 

I NCR  INCMAX  IREST  IPOST 


Input  Variables  : 


TIME 

TMAX 

TREST 

TPOST 

DTMIN 
DTMAX 
I  NCR 
INCMAX 
IREST 

IPOST 


Time  at  start  of  solution 
Maximum  time  value 

Time  interval  at  which  restart  files  are  to  be 
written  (if  IREST  >  0) 

Time  interval  at  which  results  files  are  to  be 
written  (if  IPOST  >  0) 

Minimum  time  step  value 

Maximum  time  step  value 

Increment  number  at  start  of  solution  (usually  0) 
Maximum  number  of  time  increments  to  be  performed 

Restart  save  flag  (  >0  if  restart  files  are  to  be 

crea  ted ) 

Results  save  flag  (  >0  if  results  files  are  to  be 
created  for  plotting) 


Notes  : 


(1)  -  Presently,  only  the  parameters  TIME,  TMAX,  INCR,  and 
INCMAX  are  used  when  a  coupled  analysis  is  performed 


TOLERANCES  Input  Block 


< 


I  ' 


Header  :  TOLE 

End  :  (1  Record) 

Record  Format  #1  : 

•  •  «  .  1 0  •  •  «  *20  .  .  «  *30  .  •  «  *40 

DTFRAC  V0LT0L  SDTTOL  HGDAMP 

Input  Variables  : 

DTFRAC  «  Fraction  of  the  maximum  stable  time  step  to  be 
used  in  the  solution  (  Default  -  0.75  ) 

V0LT0L  =  Relative  volume  change  to  be  permitted  before  a 
mesh  rezone  is  performed  (  Default  -  0.10  ) 

SDTTOL  =  Relative  change  in  estimated  stable  time  step  to 
be  permitted  before  a  mesh  rezone  is  performed 
(  Default  =  0.80) 

HGDAMP  =  Hourglass  damping  fraction  (  Default  -  0.3  ) 

Notes  : 

(1)  -  Although  automatic  time-stepping  is  used  in  the 
numerical  solution,  high  velocity  problems  often 
require  a  finer  step  than  the  computed  stability 
limit.  A  small  DTFRAC  will  be  required  when  the 
velocities  present  in  the  problem  approach  local 
wave  speeds  in  the  material. 
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(2)  -  Rezoning  will  be  performed  if  EITHER  the  volume 

change  or  time  step  change  tolerance  is  exceeded. 
Rezoning  is  minimized  when  VOLT0L=1.0  and 
SDTTOL-O.O.  As  VOLTOL  approaches  0  and/or  DSTTOL 
approaches  1 ,  rezoning  will  be  performed  more 
frequently  (i.e.,  for  less  extreme  distortions). 

(3)  -  By  default,  no  hourglass  damping  is  used  in  the 

solution.  When  hourglass  damping  is  used,  values 
of  HGDAMP  in  the  range  0.2  -  0.3  are  typical. 

Note  that  the  use  of  hourglass  damping  decreases 
the  allowable  time  step  by  a  mutiplicative  factor 
of  SQRT( 1 . -HGDAMP**2) -HGDAMP ;  in  fact,  for  values 
of  HGDAMP  above  0.7071 ,  the  allowable  step  size 
becomes  zero. 
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PROPERTIES  Input  Block 


Header  :  PROP 
End  :  (1  Record) 

Record  Format  #1  : 


RHO 

G 

KL 

KQ 

SY 

ET 

* 


Input  Variables  : 


RHO 

G 

KL 

KQ 

SY 

ET 


=  Reference  value  of  material  density  (pressure=0) 
-  Elastic  shear  modulus 
=  Linearized  bulk  modulus 
=  Quadratic  pressure  coefficient 
=  Yield  stress 

=  Secondary  (plastic)  tangent  modulus 


Notes  : 


(1)  -  The  material  is  assumed  to  obey  an  isotropic,  hypoelastic 

constitutive  law,  in  which  the  co-rotational  rate  of  the 

y 

Cauchy  stress,  denoted  by  o,  is  related  to  the  rate  of 
deformation  d  by 


l  -  (K  -  2G/3)  I  tr(d)  +  2G  d 
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v 

The  Jaumann  stress  rate  o  is  related  to  the  material  time 
derivative  of  o  through  the  definition 

y 

o  =  do/dt  wo  +  ow 

in  which  w  is  the  spin  tensor.  The  instantaneous  bulk 
modulus  K  is  determined  as  described  in  Note  (2)  below. 

(2)  The  equation  of  state  governing  pressure-volume  behavior 
of  the  material  is 

P  •  K  In  ( P / PQ ) 

where  p  is  the  hydrostatic  pressure  (positive  in 
compression).  The  instantaneous  bulk  modulus  is  defined 
by 

K  =  K,  (tension,  p/p  <  1) 

Li  O 

2 

K  =  K,  +  (— ^-1 )  Kn  (compression,  p/p^  >  1 ) 
u  PQ  w  o 

Therefore,  the  bulk  behavior  in  compression  exhibits  a 
"stiffening"  effect  which  does  not  occur  in  tension. 

(3)  When  the  effective  stress  in  an  element  reaches  the  yield 
value,  the  constitutive  relationship  used  is  of  the 
elastic-plastic-hydrodynamic  type.  Thus,  material  which 
deforms  permanently  at  very  low  stress  can  be  simulated 
by  specifying  a  low  value  of  SY. 
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REFERENCE  STATIONS  Input  Block 


Header  :  REFE 
End  : 

Record  Format  #1  : 

.  .  5  .  .10  .  .15 

NX  NY  NZ 

Record  Format  #2  :  (repeat  as  needed  to  define  NX  reference 

stations,  then  NY,  then  NZ) 

.  .  .  .10  .  .  .  .20  .  .  .30  .  .  .40 . 80 

REF( 1 )  REF( 2 )  REF( 3 )  REF(H)  ...  REF(8) 

REF( 9 )  REF(10)  REF(II)  REF(12)  ...  (etc.) 


Input  Variables  : 

NX/Y/Z  =  Number  of  rezone  reference  stations  in  the  X,  Y,  and 
Z  directions,  respectively 

REF(i)  =  i-th  reference  station  in  a  particular  direction.  A 

new  input  line  should  be  started  for  each  of  the  X,  Y, 
and  Z  directions.  Values  are  entered  8  per  line,  in 
ascending  order,  on  as  many  lines  as  necessary. 

Notes : 

(1)  -  Presently,  rezoning  algorithm  II  (subsection  2.6.3)  is 

used  exclusively  for  rezoning.  With  this  algorithm,  the 
only  effect  of  the  reference  station  input  is  to  establish 
constant  plotting  limits  for  the  SBI  solution,  using  the 
lowest  and  highest  coordinate  values  in  each  direction. 
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Therefore,  NX,  NY,  and  NZ  should  each  normally  be  set  to 
2.  Specify  X,  Y,  Z  reference  station  values  (2  each) 
which  define  the  limits  for  plotting  throughout  the 
solution. 


> 


* 
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RIGID 


WALL 


B 


C 


Input  Block 


Header  :  RIGI 

End  :  Blank  line 

Record  Format  #1  :  (repeat  as  necessary  to  define  all  rigid 

wall  b.c.'s  for  the  model) 

•  •  *10  .  •  . 20  .  «  *30  •  •  .40  •  •  .50  .  .  . 60  .  .  .70  .  .  .60 

ABCDEFGH 

Input  Variables  : 

A,B,...  =  Coefficients  for  a  single  rigid-wall  boundary  con¬ 
dition.  The  constraint  defined  by  coefficients  A, 

B,  ...  ,  H  has  the  form: 

Ax  +  By  +  Cz  -  (D  +  EsinFt  +  GcosHt)  >  0 

For  example,  to  define  the  condition  x>2/  let  A=l/ 
B=C=0 /  D=2 /  and  E=F=G=H=0.  To  define  a  constraint 
y < 7 ,  let  A=C=E=F=G=H=0 /  B=-l  and  D=-7.  The  E  ,  F  /G  / 
and  H  terms  permit  the  specification  of  time-vary¬ 
ing  conditions/  such  as  x  >  2  sin  4t  (A=l/  E=2  and 
F=4 )  . 


* 
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SYMMETRY 


B 


C  .  S  Input  Block 


Header  :  SYMM 

End  :  Blank  line 

Record  Format  #1  :  > 


.  .  .  5  . 

.  .10  .  .  .15 

.  .  .20  . 

.  .25 

.  .  .  30 

IDIR 

VALUE 

IFLAGI 

VELOCI 

.  .  .35  . 

.  .40  .  .  .45 

.  .  .50  . 

.  .55 

.  .  .  60 

IFLAG2 

VELOC2 

IFLAG3 

VELOC3 

Input  Variables  : 

IDIR  =  Orientation  of  plane  on  which  boundary  condition  is 

to  be  imposed  (e.g./  IDIR=1  for  an  X=constant  plane). 

VALUE  =  Position  of  plane  (value  of  X(IDIR)). 

IFLAGi  =  Switch  for  direction  "i"  for  this  boundary  condition. 

If  IFLAGi  =  0/  no  constraints  are  imposed  on  the  i-th 
component  of  velocity.  When  IFLAGi  >  0,  the  velocity 
in  the  i-th  coordinate  direction  is  VELOCi. 

VELOCi  =  Prescribed  velocity  for  coordinate  direction  "i".  The 
value  of  VELOCi  is  ignored  if  IFLAGi  =  0. 

Notes  : 

(1)  -  The  boundary  condition  defined  by  a  typical  input  line 

could  be  stated  as:  "on  the  plane  X(IDIR)  =  VALUE,  for 
any  direction  i=l,2,3  for  which  IFLAGi  >  0,  the  veloc¬ 
ity  is  equal  to  VELOCi". 

(2)  -  This  constraint  type  may  be  used  to  define  boundary 

conditions  other  than  symmetry  (although  this  is  its 
most  common  use). 
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TITLE  Input  Block 


Header  :  TITL 
End  :  (1  Record) 

Record  Format  #1  : 

•  •  .10  •  •  .20  .  .  .30  .  •  .40  .  •  ...  .  .  .30 

TITLE  (  up  to  80  characters  ) 


Input  Variables  : 

TITLE  =  Any  alphanumeric  string,  up  to  80  characters  long 
to  be  used  as  a  page  header  in  the  problem  output 


5.4  MODELING  GUIDELINES 


Finite  element  modeling  of  a  fully  coupled  soft  body  impact 
problem  is  no  more  difficult  than  modeling  of  more  conventional 
structural  dynamics  problems.  However,  the  radically  different 
methods  of  solution  employed  in  the  two  meshes  and  the  techniques 
used  to  establish  the  impact  interface  require  that  certain 
guidelines  be  followed  if  accurate  results  are  to  be  obtained. 
This  section  is  devoted  to  the  discussion  of  these  specialized 
modeling  requirements. 

First,  we  wish  to  re-emphasize  that  the  finite  elements  and 
the  solution  methodology  used  in  the  two  finite  element  meshes 
are  quite  different.  In  the  structure  mesh,  the  finite  elements 
are  usually  of  higher  order  (quadratic) ,  so  that  fewer  elements 
are  needed  for  high  accuracy.  Moreover,  the  structural  response 
is  computed  using  implicit  methods  which  involve  large  matrix 
manipulations.  This  means  that  nodal  numbering  has  an  important 
effect  upon  the  efficiency  of  the  solution,  and  measures  should 
be  taken  during  model  preparation  to  control  the  matrix 
bandwidth.  The  soft  body  mesh  is  composed  of  low-order  elements, 
and  larger  numbers  of  elements  are  typical  for  realistic 
problems;  however,  the  solution  in  this  mesh  is  explicit  (using 
small  time  steps  and  no  matrix  manipulations),  and  numbering  of 
the  mesh  has  absolutely  no  effect  upon  solution  efficiency. 

Boundary  conditions  in  the  soft  body  mesh  are  a  potential 
source  of  solution  error,  since  they  are  applied  in  a  completely 
different  manner  from  those  in  the  structural  mesh.  Situations 
in  which  extreme  care  is  necessary  usually  arise  when  symmetry  is 
used  to  reduce  the  size  of  a  problem.  When  symmetry  conditions 
must  be  imposed  in  the  soft  body  mesh,  especially  in  the  vicinity 
of  an  impact  site,  the  following  procedures  are  suggested: 

•  apply  BOTH  symmetry  ( SYMM )  -and  rigid-wall  (RIGI)  boundary 
conditions  along  the  symmetry  edges  of  the  soft-body 
mesh;  and 
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•  create  additional  interface  elements  which  cross  the 
symmetry  boundaries  in  the  vicinity  of  an  impact. 

The  first  measure  is  desirable  because  boundary  conditions 
in  the  soft  body  mesh  are  applied  not  to  specific  nodes,  but  to 
all  nodes  which  instantaneously  lie  on  a  specified  surface.  The 
RIGI  (rigid— wall)  conditions  prevent  nodes  from  crossing  the 
symmetry  plane,  a  situation  which  may  occur  when  imposing  the 

c  interface  conditions  in  a  contact  zone;  the  SYMM  (symmetry) 

conditions  impose  the  desired  zero-velocity  conditions.  If  this 

*  procedure  is  not  followed,  the  following  sequence  of  events  is 
possible : 

(a)  a  node  "penetrates"  the  structure  mesh,  and  must  be 

moved  in  order  to  enforce  compatibility  of  the  two 
meshes ; 

(b)  when  the  node  is  moved  (always  normal  to  the  surface 

of  the  structure  mesh),  it  is  pushed  off  the  plane  of 
symmetry,  usually  to  the  side  opposite  the  model; 

(c)  the  desired  symmetry  conditions  are  not  imposed,  since 

the  node  no  longer  lies  on  the  symmetry  plane;  and 

(d)  in  some  cases,  the  subsequent  motion  of  the  node  is 

such  that  it  may  actually  pass  the  edge  of  the 
structure  mesh. 

The  addition  of  a  RIGId  wall  condition  causes  the  node  in 
question  to  be  moved  back  to  the  symmetry  plane  immediately 
following  step  (b)  above,  and  both  the  symmetry  and  contact 
boundary  conditions  will  be  satisfied. 

,  The  second  modeling  device  mentioned  above  in  connection 

with  symmetry  surface  modeling  involves  the  addition  of  interface 

•  elements  which  actually  cross  the  symmetry  plane.  This  measure 
ensures  that  the  phenomenon  mentioned  in  (d)  above  does  not  occur 
in  pathological  situations.  The  tolerances  used  in  determining 
contact  are  necessarily  very  small,  since  making  them  larger  will 
result  in  certain  nodes  applying  the  same  forces  twice  to  two  or 
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more  neighboring  elements.  In  some  instances,  these  small 
tolerances  may  cause  true  contact  conditions  on  symmetry  edges  to 
be  skipped.  The  suggested  procedure  is  to  create  an  additional 
"layer"  of  nodes  in  the  structure  mesh,  with  their  motion  tied  to 
the  motion  of  the  original  structural  nodes  by  linear  constraint 
equations.  The  extra  layer  of  nodes  then  allows  an  additional 
boundary  "layer"  of  interface  elements  to  be  specified,  through 
which  the  soft-body  mesh  may  not  pass. 

A  final  consideration  involves  controlling  the  solution  time 
step  in  the  soft  body  mesh  in  moderate  to  high  velocity  problems. 
Normally,  the  time  step  in  this  mesh  is  selected  automatically, 
based  on  estimates  of  the  maximum  wave  speed  in  the  body.  When 
impact  conditions  are  specified,  however,  nodes  in  the  soft  body 
mesh  must  be  adjusted  to  eliminate  any  "penetration"  of  the 
structure  mesh,  with  velocities  and  accelerations  on  the  impact 
surface  being  adjusted  accordingly.  When  the  product  of  the 
relative  velocity  in  the  two  meshes  and  the  time  step  in  the  soft 
body  mesh  is  larger  than  the  soft  body  mesh  spacing,  these  nodal 
coordinate  corrections  may  produce  elements  with  negative  volume 
or,  at  least,  severe  distortion  within  a  single  time  step.  In 
such  situations,  it  is  appropriate  to  specify  a  maximum  time  step 
for  the  soft  body  mesh  to  avoid  problems  of  this  type.  A 
reasonable  estimate  of  the  maximum  time  step  can  be  obtained  by 
comparing  the  initial  (usually  uniform)  velocity  of  the  soft  body 
mesh  with  the  smallest  mesh  division  in  the  same  direction;  a 
limiting  time  step  of  about  10-20%  of  the  mesh  spacing  divided  by 
the  initial  velocity  is  usually  adequate. 


4 
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5.5  ANALYSIS  RESTART  PROCEDURES  FOR  COUPLED  SOLUTION 


Full  restart  capabilities  are  available  for  the  coupled 
solution  program.  Restart  checkpoints  are  created  automatically 
whenever  a  coupled  solution  is  selected;  saving  the  restart  data 
simply  requires  that  the  restart  data  file  be  saved  on  disk  when 
a  segment  of  the  analysis  is  complete.  This  "new"  restart  file 
is  called  "HYDNEW"  in  all  machine  versions  of  the  code;  default 
file  extensions  for  FORTRAN  data  files  are  used  on  machines  such 
as  the  VAX- 11/7 80. 

Restart  of  an  interrupted  solution  is  requested  whenever 
the  MAGNA  "RESTART"  option  is  invoked.  Restart  data  is  read  from 
the  existing  file  "HYDOLD"  (as  with  HYDNEW,  default  file 
extensions  are  used  where  applicable) .  On  CDC  and  CRAY 
computers,  the  HYDOLD  file  should  be  copied  from  permanent  file 
storage  to  a  local  file  in  the  input  job  stream. 

Occasionally,  it  may  be  convenient  to  restart  a  coupled 
solution  as  a  normal  nonlinear  dynamic  analysis,  omitting  the 
hydrodynamic  mesh  entirely.  An  example  is  a  hydrodynamic  impact 
case  in  which  the  impact  loading  phase  is  complete,  and  the  two 
meshes  have  separated;  analysis  of  the  post-impact  vibration  of 
the  structural  model  does  not  require  the  use  of  the  coupled 
solution  option.  To  accomplish  such  a  restart,  simply  restart 
the  structural  analysis  as  a  normal  nonlinear  dynamic  problem 
(i.e.,  changing  IOPT(18)  to  zero).  The  files  FSDATA  and  HYDOLD 
need  not  be  supplied  as  input. 

Note  that  the  opposite  case  (restarting  a  normal  analysis 
as  a  coupled  one)  is  presently  not  permitted.  Once  a  solution 
(or  a  segment  of  a  solution)  has  been  performed  without  the 
coupling  option,  it  may  be  restarted  only  as  an  uncoupled 
analysis . 


99 


SECTION  6 

SUMMARY  AND  CONCLUSIONS 

The  soft  body  impact  analysis  described  in  this  report 
represents  an  attempt  to  treat  this  important  class  of  problems 
without  adopting  heuristic  or  highly  idealized  models  of  impact 
loading.  As  such,  it  will  find  application  in  situations  where 
impact  loading  and  the  resulting  response  are  highly  coupled  and 
cannot  be  determined  independently.  Examples  include  bird 
impacts  on  aircraft  transparencies  and  foreign  object  ingestion 
by  air-breathing  propulsion  systems. 

The  analytical  framework  developed  in  this  report  is  general 
in  scope,  and  can  be  enhanced  as  experience  in  analyzing 
realistic  problems  accumulates.  Thus  far,  our  experience  in 
simulating  real  soft-body  impacts  is  limited,  and  correlation  of 
analytical  and  experimental  results  has  been  impeded  by  a  lack  of 
mechanical  properties  data.  However,  several  of  the  modeling 
exercises  performed  thus  far  have  provided  useful  knowledge  of  a 
pragmatic  nature.  In  particular,  the  following  observations  can 
be  made  based  upon  our  present  experience: 

o  with  a  fully  three-dimensional  model  based  on  "primitive” 
variables,  a  slight  amount  of  shear  stiffness  (or  fluid 
viscosity)  is  always  necessary  for  a  stable  solution  which 
is  free  from  spurious  modes  of  deformation; 

o  the  ant i-hourglass ing  technique  suggested  in  Reference 
[21]  is  quite  reliable,  and  moderate  hourglass  damping 
ratios  (0.2-0. 3)  are  virtually  always  adequate; 

o  an  elastic-plastic-failure  model  is  somewhat  easier  to 
work  with  in  soft-body  impact  problems  than  a  Newtonian 
fluid  model,  since  the  fluid  model  parameters  are 
difficult  to  obtain  and  since  the  physical  situation 
normally  involves  at  least  some  elastic  behavior; 
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o  the  use  of  a  nonlinear  model  for  bulk  behavior  of  the  soft 
body  material  can  produce  predicted  results  which  are 
radically  different  (and  usually  improved)  from  a  linear 
pressure-volume  model; 

o  the  application  of  symmetry  boundary  conditions  in 
combination  with  impact  constraints  requires  careful 
modeling  to  maintain  correct  behavior  of  the  hydrodynamic 
mesh  (subsection  5.4); 

o  time  steps  in  the  soft  body  mesh  can  be  computed  quite 
accurately  from  the  instantaneous  properties  of  individual 
finite  elements,  permitting  continuous  and  automatic 
adjustment  of  the  solution  step  size  (subsection  2.4); 

o  modification  of  the  computed  time  step  for  the  soft  body 
mesh  normally  can  be  limited  to  specifying  a  maximum  time 
step,  when  the  impact  velocity  is  large  enough  that  some 
elements  might  be  compressed  by  their  entire  length  within 
a  single  time  step  upon  impact  (subsection  5.4); 

o  force  information  obtained  from  an  impact  solution  may 
become  inaccurate  when  too  few  steps  are  used  in  the 
structural  mesh,  since  peak  forces  which  occur  early  in 
the  impact  event  may  build  and  decay  within  a  single  time 
step  (Section  3) • 


The  next  important  step  toward  the  routine  analysis  of 
practical  problems  in  soft  body  impact  is  the  validation  of  the 
present  methodology  using  experimental  data  from  impact  tests  on 
deformable  structures.  Full  scale  impact  test  results  for 
transparency  birdstrikes  exist  [22],  and  can  be  used  for  the 
validation.  The  collection  of  mechanical  properties  data  for 
both  real  and  substitute  bird  materials  in  support  of  this 
validation  activity  is  recommended  as  well. 
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